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For want of a nail the shoe was lost. 
For want of a shoe the horse was lost. 
For want of a horse the rider was lost. 
For want of a rider the battle was lost. 
For want of a battle the kingdom was lost. 

And all for the want of a horseshoe nail. 
For Want of a Nail (proverbial rhyme) 



Summary. We present a survey of the theory of the Lyapunov Characteristic Expo- 
nents (LCEs) for dynamical systems, as well as of the numerical techniques developed 
for the computation of the maximal, of few and of all of them. After some histor- 
ical notes on the first attempts for the numerical evaluation of LCEs, we discuss 
in detail the multiplicative ergodic theorem of Oseledec [102| , which provides the 
theoretical basis for the computation of the LCEs. Then, we analyze the algorithm 
for the computation of the maximal LCE, whose value has been extensively used as 
an indicator of chaos, and the algorithm of the so-called 'standard method', devel- 
oped by Benettin et al. [I4], for the computation of many LCEs. We also consider 
different discrete and continuous methods for computing the LCEs based on the QR 
or the singular value decomposition techniques. Although, we are mainly interested 
in finite-dimensional conservative systems, i. e. autonomous Hamiltonian systems 
and symplectic maps, we also briefiy refer to the evaluation of LCEs of dissipative 
systems and time series. The relation of two chaos detection techniques, namely the 
fast Lyapunov indicator (FLI) and the generalized alignment index (GALI), to the 
computation of the LCEs is also discussed. 

Key words: Lyapunov exponents; Multiplicative ergodic theorem; Numerical tech- 
niques; Dynamical systems; Chaos; Variational equations; Tangent map; Chaos de- 
tection methods 
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1 Introduction 



One of the basic information in understanding the behavior of a dynamical 
system is the knowledge of the spectrum of its Lyapunov Characteristic Expo- 
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nents (LCEs). The LCEs are asymptotic measures characterizing the average 
rate of growth (or shrinking) of small perturbations to the solutions of a dy- 
namical system. Their concept was introduced by Lyapunov when studying 
the stability of nonstationary solutions of ordinary differential equations [96j , 
and has been widely employed in studying dynamical systems since then. The 
value of the maximal LCE (mLCE) is an indicator of the chaotic or regu- 
lar nature of orbits, while the whole spectrum of LCEs is related to entropy 
(Kolmogorov-Sinai entropy) and dimension-like (Lyapunov dimension) quan- 
tities that characterize the underlying dynamics. 

By dynamical system we refer to a physical and/or mathematical system 
consisting of a) a set of I real state variables xi,X2 ...jXi, whose current 
values define the state of the system, and b) a well-defined rule from which 
the evolution of the state with respect to an independent real variable (which 
is usually referred as the time t) can be derived. We refer to the number / of 
state variables as the dimension of the system, and denote a state using the 
vector X = {xi,X2 ■ ■ ■ , Xi), or the matrix x = [xi Xi . . . xi\^ notation, where 
C^) denotes the transpose matrix. A particular state x corresponds to a point 
in an /-dimensional space 5, the so-called phase space of the system, while a 
set of states x(i) parameterized by t is referred as an orbit of the dynamical 
system. 

Dynamical systems come in essentially two types: 

1. Continuous dynamical systems described by differential equations of the 
form 

with dot denoting derivative with respect to a continuous time t and f 
being a set of / functions /i , /2 ■ ■ ■ , fi known as the vector field. 

2. Discrete dynamical systems or maps, described by difference equations of 
the form 

x„+i = f(x„), 

with f being a set of I functions /i , /2 ■ ■ ■ , fi and x„ denoting the vector 
X at a discrete time t = n (integer). 

Let us now define the term chaos. In the literature there are many defini- 
tions. A brief and concise presentation of them can be found for example in 
[90] . We adopt here one of the most famous definitions of chaos due to De- 
vaney [551 P- 50], which is based on the topological approach of the problem. 



Definition 1. Let V be a set and f : V V a map on this set. We say that 
f is chaotic on V if 

1. f has sensitive dependence on initial conditions. 

2. f is topologically transitive. 

3. periodic points are dense in V . 
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Let us explain in more detail the hypothesis of this definition. 

Definition 2. f . V ^ V has sensitive dependence on initial conditions 

if there exists 6 > such that, for any x G and any neighborhood A of x, 
there exist y G A and n > 0, such that |f"(x) — f"(y)| > S, where f" denotes 
n successive applications off. 

Practically this definition implies that there exist points arbitrarily close to 
X which eventually separate from x by at least 6 under iterations of f . We 
point out that not all points near x need eventually move away from x under 
iteration, but there must be at least one such point in every neighborhood of 

X. 

Definition 3. f : V ^ V is said to be topologically transitive if for any 

pair of open sets U, W C V there exists n > such that i"" (U) H W ^ 9 . 

This definitions implies the existence of points which eventually move under 
iteration from one arbitrarily small neighborhood to any other. Consequently, 
the dynamical system cannot be decomposed into two disjoint invariant open 
sets. 

From Definition[l]we see that a chaotic system possesses three ingredients: 

a) unpredictability because of the sensitive dependence on initial conditions, 

b) indecomposability, because it cannot be decomposed into noninteracting 
subsystems due to topological transitivity, and c) an element of regularity 
because it has periodic points which are dense. 

Usually, in physics and applied sciences, people focus on the first hypoth- 
esis of Definition [T] and use the notion of chaos in relation to the sensitive 
dependence on initial conditions. The most commonly employed method for 
distinguishing between regular and chaotic motion, which quantifies the sen- 
sitive dependence on initial conditions, is the evaluation of the mLCE xi- If 
Xi > the orbit is chaotic. This method was initially developed at the late 
70's based on theoretical results obtained at the end of the 60's. 

The concept of the LCEs has been widely presented in the literature from a 
practical point of view, i. e. the description of particular numerical algorithms 
for their computation [54l |44l [62l [92l [36] . Of course, there also exist theoret- 
ical studies on the LCEs, which are mainly focused on the problem of their 
existence, starting with the pioneer work of Oseledec |102j . In that paper the 
Multiplicative Ergodic Theorem (MET), which provided the theoretical basis 
for the numerical computation of the LCEs, was stated and proved. The MET 
was the subject of several theoretical studies afterwards |108[|114l [7^1141] . A 
combination of important theoretical and numerical results on LCEs can be 
found in the seminal papers of Benettin et al. |131 114j , written almost 30 years 
ago, where an explicit method for the computation of all LCEs was developed. 

In the present report we focus our attention both on the theoretical frame- 
work of the LCEs, as well as on the numerical techniques developed for their 
computation. Our goal is to provide a survey of the basic results on these 
issues obtained over the last 40 years, after the work of Oseledec |102] . To 
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this end, we present in detail the mathematical theory of the LCEs and dis- 
cuss its significance without going through tedious mathematical proofs. In 
our approach, we prefer to present the definitions of various quantities and to 
state the basic theorems that guarantee the existence of the LCE, citing at 
the same time the papers where all the related mathematical proofs can be 
found. We also describe in detail the various numerical techniques developed 
for the evaluation of the maximal, of few or even of all LCEs, and explain 
their practical implementation. We do not restrict our presentation to the so- 
called standard method developed by Benettin et al. [H], as it is usually done 
in the literature (see e. g. [SH |44l |92] ) , but we include in our study modern 
techniques for the computation of the LCEs like the discrete and continu- 
ous methods based on the singular value decomposition (SVD) and the QR 
decomposition procedures. 

In our analysis we deal with finite-dimensional dynamical systems and 
in particular with autonomous Hamiltonian systems and symplectic maps 
defined on a compact manifold, meaning that we exclude cases with escapes in 
which the motion can go to infinity. We do not consider the rather exceptional 
cases of completely chaotic systems and of integrable ones, i. e. systems that 
can be solved explicitly to give their variables as single-valued functions of 
time, but we consider the most general case of 'systems with divided phase 
space' [30l p. 19] for which regulail^ (quasiperiodic) and chaotic orbits co-exist. 
In such systems one sees both regular and chaotic domains. But the regular 
domains contain a dense set of unstable periodic orbits, which are followed by 
small chaotic regions. On the other hand, the chaotic domains contain stable 
periodic orbits that are followed by small islands of stability. Thus, the regular 
and chaotic domains are intricately mixed. However, there are regions where 
order is predominant, and other regions where chaos is predominant. 

Although in our report the theory of LCEs and the numerical techniques 
for their evaluation are presented mainly for conservative systems, i.e. system 
that preserve the phase space volume, these techniques are not valid only for 
such models. For completeness sake, we also briefly discuss at the end of the 
report the computation of LCEs for dissipative systems, for which the phase 
space volume decreases on average, and for time series. 

We tried to make the paper self-consistent by including definitions of the 
used terminology and brief overviews of all the necessary mathematical no- 
tions. In addition, whenever it was considered necessary, some illustrative 
examples have been added to the text in order to clarify the practical imple- 
mentation of the presented material. Our aim has been to make this review 
of use for both the novice and the more experienced practitioner interested 
in LCEs. To this end, the reader who is interested in reading up on detailed 
technicalities is provided with numerous signposts to the relevant literature. 

Throughout the text bold lowercase letters denote vectors, while matri- 
ces are represented, in general, by capital bold letters. We also note that the 



^ Regular orbits are often called ordered orbits (see e. g. [301 p. 18]). 
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most frequently used abbreviations in the text are: LCE(s), Lyapunov Char- 
acteristic Exponent(s); p-LCE, Lyapunov Characteristic Exponent of order 
p; mLCE, maximal Lyapunov Characteristic Exponent; p-mLCE, maximal 
Lyapunov Characteristic Exponent of order p; MET, multiplicative ergodic 
theorem; SVD, Singular Value Decomposition; PSS, Poincare surface of sec- 
tion; FLI, fast Lyapunov indicator; GALI, generalized alignment index. 
The paper is organized as follows: 

In Section [2] we present the basic concepts of Hamiltonian systems and 
symplectic maps, emphasizing on the evolution of orbits, as well as of deviation 
vectors about them. In particular, we define the so-called variational equations 
for Hamiltonian systems and the tangent map for symplectic maps, which 
govern the time evolution of deviation vectors. We also provide some simple 
examples of dynamical systems and derive the corresponding set of variational 
equations and the corresponding tangent map. 

Section[3] contains some historical notes on the first attempts for the appli- 
cation of the theoretical results of Oseledec |102j for the actual computation 
of the LCEs. We recall how the notion of exponential divergence of nearby 
orbits was eventually quantified by the computation of the mLCE, and we 
refer to the papers where the mLCE or the spectrum of LCEs were computed 
for the first time. 

The basic theoretical results on the LCEs are presented in Section|4]follow- 
ing mainly the milestone papers of Oseledec |102j and Benettin et al. [T3lll4j. 
In Section 14.11 the basic definitions and theoretical results of LCEs of various 
orders are presented. The practical consequences of these results on the com- 
putation of the LCEs of order 1 and of order p > 1 are discussed in Sections 
[421 and gj] respectively. Then, in Section gj] the MET of Oseledec [102] is 
stated in its various forms, while its consequences on the spectrum of LCEs 
for conservative dynamical systems are discussed in Section [4.51 

Section [5] is devoted to the computation of the mLCE xij which is the 
oldest chaos indicator used in the literature. In Section [STT] the method for the 
computation of the mLCE is discussed in great detail and the theoretical basis 
of its evaluation is explained. The corresponding algorithm is presented in 
Section [5. 2[ while the behavior of xi for regular and chaotic orbits is analyzed 
in Section [5T3] 

In Section [6] the various methods for the computation of part or of the 
whole spectrum of LCEs are presented. In particular, in Section (HH] the stan- 
dard method developed in [119[ [T4] . is presented in great detail, while the 
corresponding algorithm is given in Section 16.21 In Section 16.31 the connec- 
tion of the standard method with the discrete QR decomposition technique 
is discussed and the corresponding QR algorithm is given, while Section 16.41 
is devoted to the presentation of other techniques for computing few or all 
LCEs, which are based on the SVD and QR decomposition algorithms. 

In Section [7] we briefly refer to various chaos detection techniques based 
on the analysis of deviation vectors, as well as to a second category of chaos 
indicators based on the analysis of the time series constructed by the coordi- 
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nates of the orbit under consideration. The relation of two chaos indicators, 
namely the fast Lyapunov indicator (FLI) and the generalized alignment index 
(GALI), to the computation of the LCEs is also discussed. 

Although the main topic of our presentation is the theory and the com- 
putation of the LCEs for conservative dynamical systems, in the last section 
of our report some complementary issues related to other types of dynamical 
systems are concisely presented. In particular, Section 18.11 is devoted to the 
computation of the LCEs for dissipative systems, while in Section 18.21 some 
basic features on the numerical computation of the LCEs from a time series 
are presented. 

Finally, in the appendix |^ we present some basic elements of the exterior 
algebra theory in connection to the evaluation of wedge products, which are 
needed for the computation of the volume elements appearing in the defini- 
tions of the various LCEs. 



2 Autonomous Hamiltonian systems and symplectic 
maps 

In our study we consider two main types of conservative dynamical systems: 

1. Continuous systems corresponding to an autonomous Hamiltonian system 
of N degrees {NT)) of freedom having a Hamiltonian function 



H{qi,q2, ■ ■ ■,qNiPi,P2, ■ ■ ■ ,Pn) = h = constant, 



(1) 



where Qi and Pi, i = 1,2, . . . , N are the generalized coordinates and con- 
jugate momenta respectively. An orbit in the I — 2A^-dimensional phase 
space S of this system is defined by a vector 

x(0 = {qi{t),q2{t), ■ ■ . ,qN{t),Pi{t),P2{t), ■ ■ ■ ,PN{t)), 

with Xi = qi, Xij^js! — Pi, i — 1,2, . . . , N . The time evolution of this orbit 
is governed by the Hamilton equations of motion, which in matrix form 
are given by 



X = f (x) = 



OH 
dp 



OH 

^c^ 



= •i2N 



DH, 



(2) 



with q = {qi{t),q2{t), . . . ,gAr(t)), p = {pi{t),p2(t), . . . ,PN{t)), and 



DH 



dH dH 
dqi dq2 



dH OH dH 
dqN dpi dp2 



dH 
dpN 



Matrix 32N has the following block form 

32N = 



On In 
-In On 
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with Ijv being the N x N identity matrix and Ojv being the N x N matrix 
with all its elements equal to zero. The solution of Q is formally written 
with respect to the induced flow ^* : — > 5 as 

x(i)=<P*(x(0)). (3) 

2. Symplcctic maps of I = 2N dimensions having the form 

x„+i = f(x„). (4) 

A symplectic map is an area-preserving map whose Jacobian matrix 

Of: 



M = Df(x) = |^ 



dfi dfi 

dxi dx2 

dh, dh 

dxi dx2 



df2N df2N 

dxi dx2 



dx 



dx2 



dhN 

dx2N 



satisfies 

M"^ • Jaw ■ M = 32N- 

The state of the system at the discrete time t = n is given by 

x„ = (xo) = (f)" (xo) , 
where (f)" (xo) ~ f (f (• • • f (xq) • • •)), n times. 



(5) 
(6) 



2.1 Variational equations and tangent map 

Let us now turn our attention to the (continuous or discrete) time evolution 
of deviation vectors w from a given reference orbit of a dynamical system. 
These vectors evolve on the tangent space T^S of S. We denote by dx?^* the 
linear mapping which maps the tangent space of S at point x onto the tangent 
space at point ??*(x), and so we have dy^'P* : ^^'(x)*^ with 

w(0 =dx<?*w(0), (7) 

where w(0), w(t) are deviation vectors with respect to the given orbit at times 
t — Q and t > respectively. 

In the case of the Hamiltonian system ([1]) an initial deviation vector 
w(0) = (5a;i(0), 52:2(0), (5a;2Ar(0)) from the solution x(t) ^ evolves on 
the tangent space 7^5 according to the so-called variational equations 

(9f 

W = Df(x(t)) ■ w = ^(x(t)) • w = [32N ■ D2H(x(i))] • w =: A{t) • w , (8) 

with D^H(x(t)) being the Hessian matrix of Hamiltonian ([T]) calculated on 
the reference orbit x(<) fS]), i. e. 
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D2H(x(t)),j - 



dxjdxi 



, t,j^l,2,...,2N. 

<Z>*(x(0)) 

We underline that equations ^ represent a set of linear differential equations 
with respect to w, having time dependent coefficients since, matrix A(i) de- 
pends on the particular reference orbit, which is a function of time t. The 
solution of ([8]) can be written as 

w(t) = Y{t) ■ w(0), (9) 

where Y(i) is the so-called fundamental matrix of solutions of ([5]) , satisfying 
the equation 

Y(t) = Df(x(i)) • Y{t) = A{t) ■ Y{t) , with Y(0) hN- (10) 

In the case of the symplectic map ^ the evolution of a deviation vector 
w„ , with respect to a reference orbit x„ , is given by the corresponding tangent 
map 

dt 

w„+i = Df(x„) • w„ = ^(x^ 



■ w„ ~. M„ • w„. 

Thus, the evolution of the initial deviation vector wq is given by 

w„ = M„_i ■ M„_2 • . . . • Mo • v^fQ =: Y„ • wo, 
with Y„ satisfying the relation 

Y„+i = M„ • Y„ = Df(x„) • Y„, with Yo = l2N- 

2.2 Simple examples of dynamical systems 



(11) 
(12) 

(13) 



As representative examples of dynamical systems we consider a) the well- 
known 2D Henon-Heiles system |72], having the Hamiltonian function 



H2 = ]:{pI+vI) + + ?/) + x'y - \y\ 



with equations of motion 



X = 



X 

y 

Px 

Py 



J4 DH, 





X + 2xy 






= Px 


J4 


y + x'^ — 




y 


^Py 


Px 




Px 


= —X 




Py 




.Py 


= -y 



2xy 



(14) 



(15) 



X 



(mod27r), (16) 



and b) the 4-dimensional (4d) symplectic map 

X2,n+1 = 2^2, n + 3^4,™ 

X3,n+i = xs^n " sin(a;i,„+i ) - /i[l - cos(a;i,„+i + X2.„+i)] 
Xi,n+i = " K sin(a;2,„+i ) - m[1 - cos(a;i^„+i + X2,„+i)] 

with parameters v, k and /u. All variables are given (mod 27r), so Xi^n G [tt, tt), 
for i — 1,2, 3, 4. This map is a variant of Froeschle's 4d symplectic map [55] 
and its behavior has been studied in [3Tlll23| . It is easily seen that its Jacobian 
matrix satisfies equation ([5|). 
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2.3 Numerical integration of variational equations 



When dealing with Hamiltonian systems the variational equations ([5]) have to 
be integrated simultaneously with the Hamilton equations of motion Let 
us clarify the issue by looking to a specific example. The variational equations 
of the 2D Hamiltonian (fT4l) are 



Sx 




10" 




Sx 


Sy 




1 




Sy 






-\-1y -2x 




Spx 


JPy. 




-2x -1 + 27/0 




_Spy_ 


f 


Sx - 


= 6px 





(17) 



Sy = Spy 
SPx - (-1 



2y)Sx + i-2x)Sy 



[Spy = i^2x)Sx + {-1 + 2y)Sy 

This system of differential equations is linear with respect to Sx^ Sy, Spx, 
Spy, but it cannot be integrated independently of system ((T5)) since the x 
and y variables appear explicitly in it. Thus, if we want to follow the time 
evolution of an initial deviation vector w(0) with respect to a reference orbit 
with initial condition x(0), we are obliged to integrate simultaneously the 
whole set of differential equations and P?|) . 

A numerical scheme for integrating the variational equations ([8]), which 
exploits their linearity and is particularly useful when we need to evolve more 
than one deviation vectors is the following. Solving the Hamilton equations of 
motion Q by any numerical integration scheme we obtain the time evolution 
of the reference orbit In practice this means that we know the values x(ii) 
for ti = i At, i = 0, 1, 2, . . ., where At is the integration time step. Inserting 
this numerically known solution to the variational equations ^ we end up 
with a linear system of differential equations with constant coefficients for 
every time interval [ti, ti + At), which can be solved explicitly. 

For example, in the particular case of Hamiltonian (fH)) . the system of 
variational equations p?)) becomes 



iorte[t^,t, + At), (18) 



Sx = Spx 
Sy = Spy 

Sp, = [-1 - 2y{t^)] Sx + [~2x{t,)] Sy 
Spy = [-2x{U)] Sx + [-1 + 2y{t,)] Sy 

which is a linear system of differential equations with constant coefficients 
and thus, easily solved. In particular, equations (fTS)) can by considered as the 
Hamilton equations of motion corresponding to the Hamiltonian function 



Hspl + spl) + Hi^ 



Hv{Sx, Sy, Spx, Spy) = 
2y{ti)]Sx'' + [l^2y{t,)] Sy'' 



2 [2x{ti)] SxSy] 



(19) 
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The Hamiltonian formalism (|19p of the variational equations (1181) is a 
specific example of a more general result. In the case of the usual Hamiltonian 
function 



1 ^ 



(20) 



with V{q) being the potential function, the variational equations ^ for the 
time interval [ti,ti + At) take the form (see e. g. [12j ') 



J AT 



N 





6q 







with = {6qi{t),Sq2{t), . . .,SqN{t)), 6p = (Spi{t),6p2{t) . 



D^V(q(<,)),fc = 



jdqk 



, j,k = 1,2,. 



,SpN{t)), and 



q(t,) 



Thus, the tangent dynamics of (|20p is represented by the Hamiltonian function 
(see e. g. [105]) 

N N 



2.4 Tangent dynamics of symplectic maps 



In the case of symplectic maps, the dynamics on the tangent space, which is 
described by the tangent map pip , cannot be considered separately from the 
phase space dynamics determined by the map ([4]) itself. This is because the 
tangent map depends explicitly on the reference orbit x„. 
For example, the tangent map of the 4d map (fT6|) is 

^a;i,ri+i 

SX2,n+l 

Sx4.n+1 

with 

an = -:/COs(xi,„+i) - ^Sin(xi^„+i + X2^n+l) 

b„ = -/isin(a;i,„+i + X2,n+i) 

Cn = -KCOs(a;2,„+i) - /isin(a::i,„+i + X2.,n+i) 

which explicitly depend on xi,„, X2^m x^^m a^4,n- Thus, the evolution of a 
deviation vector requires the simultaneous iteration of both the map (I16p and 
the tangent map ([2T|) . 



= OXiji + dXz,n 
= Sx2,n + Sx4^n 

= QnSxi^n + bnSx2.n + (1 + anjSx^^n + bnSx4. 
= bnSxi,n + CnSx2.n + bnSxz.n + (1 + Cn)5xi,, 
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3 Historical introduction: The early days of LCEs 

Prior to the discussion of the theory of the LCEs and the presentation of the 
various algorithms for their computation, it would be interesting to go back 
in time and see how the notion of LCEs, as well as the nowadays taken for 
granted techniques for evaluating them, were formed. 

The LCEs are asymptotic measures characterizing the average rate of 
growth (or shrinking) of small perturbations to the orbits of a dynamical 
system, and their concept was introduced by Lyapunov [96| . Since then they 
have been extensively used for studying dynamical systems. As it has already 
been mentioned, one of the basic features of chaos is the sensitive dependence 
on initial conditions and the LCEs provide quantitative measures of response 
sensitivity of a dynamical system to small changes in initial conditions. For a 
chaotic orbit at least one LCE is positive, implying exponential divergence of 
nearby orbits, while in the case of regular orbits all LCEs are zero. Therefore, 
the presence of positive LCEs is a signature of chaotic behavior. Usually the 
computation of only the mLCE xi is sufficient for determining the nature of 
an orbit, because xi > guarantees that the orbit is chaotic. 

Characterization of the chaoticity of an orbit in terms of the divergence 
of nearby orbits was introduced by Henon and Heiles [75] and further used 
by several authors (e. g. [H |5T1 [5l [ml [22l |2T] ) . In these studies two initial 
points were chosen very close to each other, having phase space distance of 
about 10-'^ - 10-^ and were evolved in time. If the two initial points were 
located in a region of regular motion their distance increased approximately 
linearly with time, while if they were belonging to a chaotic region the distance 
exhibited an exponential increase in time (Figure [1]) . 

Although the theory of LCEs was applied to characterize chaotic motion 
by Oseledec [102] . quite some time passed until the connection between LCEs 
and exponential divergence was made clear [10|, 1106] . It is worth mentioning 
that Casartelli et al. [H] defined a quantity, which they called 'stochastic 
parameter', in order to quantify the exponential divergence of nearby orbits, 
which was realized afterwards in [TiT to be an estimator of the mLCE for 
t —>■ oo. 

So, the mLCE xi was estimated for the first time in [10] . as the limit 
for i CX3 of an appropriate quantity Xi{t), which was obtained from the 
evolution of the phase space distance of two initially close orbits. In this paper 
some nowadays well-established properties of Xi{t) were discussed, like for 
example the fact that Xi {t) tends to zero in the case of regular orbits following 
a power law oc t~^, while it tends to nonzero values in the case of chaotic orbits 
(Figure [2]). The same algorithm was immediately applied for the computation 
of the mLCE of a dissipative system, namely the Lorenz system 

The next improvement of the computational algorithm for the evaluation 
of the mLCE was introduced in [34] , where the variational equations were used 
for the time evolution of deviation vectors instead of the previous approach 
of the simultaneous integration of two initially close orbits. This more direct 
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Fig. 1. Typical behavior of the time evolution of the distance D between two 
initially close orbits in the case of regular and chaotic orbits. The particular results 
are obtained for a 2D Hamiltonian system describing a Toda lattice of two particles 
with unequal masses (see |22) for more details). The initial Euclidian distance of the 
two orbits in the 4-dimensional phase space is Do — 10~®. D exhibits a linear (on 
the average) growth when the two orbits are initially located in a region of regular 
motion (left panel), while it grows exponentially in the case of chaotic orbits (right 
panel). The big difference in the values of D between the two cases is evident since 
the two panels have the same horizontal (time) axis but different vertical ones. In 
particular, the vertical axis is linear in the left panel and logarithmic in the right 
panel (after [22]). 

approach constituted a significant improvement for the computation of the 
mLCE since it allowed the use of larger integration steps, diminishing the real 
computational time and also eliminated the problem of choosing a suitable 
initial distance between the nearby orbits. 

In [TT] a theorem was formulated, which led directly to the development 
of a numerical technique for the computation of some or even of all LCEs, 
based on the time evolution of more than one deviation vectors, which are kept 
linearly independent through a Gram-Schmidt orthonormalization procedure 
(see also [5]). This method was explained in more detail in [119j . where it 
was applied to the study of the Lorenz system and was also presented in [12j . 
where it was applied to the study of an A^D Hamiltonian system with A'' 
varying from 2 to 10. 

The theoretical framework, as well as the numerical method for the com- 
putation of the maximal, some or even all LCEs were given in the seminal 
papers of Benettin et al. [131 [13] . In [T3] the complete set of LCEs was cal- 
culated for several different Hamiltonian systems, including four and six di- 
mensional maps. In Figure [3] we show the results of [14] concerning the 3D 
Hamiltonian system of [33] ■ The importance of the papers of Benettin et al. 
[T31 [Tl] is reflected by the fact that almost all methods for the computation of 
the LCEs are more or less based on them. Immediately the ideas presented in 
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Fig. 2. Evolution of Xi{t) (denoted as fc„) with respect to time t (denoted by 
n X r) in log-log scale for several orbits of the Henon-Heiles system (I14p . In the 
left panel Xi{t) is computed for 5 different regular orbits at different energies H2 
(denoted as E) and it tends to zero following a power law oc . A dashed straight 
hne corresponding to a function proportional to is also plotted. In the right 
panel the evolution of Xi (t) is plotted for three regular orbits (curves 1-3) and 
three chaotic ones (curves 4-6) for H2 = 0.125. Note that the values of the initial 
conditions given in the two panels correspond to gi — q^. — y, pi — Px, P2 = Py 
in dill (after [TO]). 

p!3l [T4] were used for the computation of the LCEs for a variety of dynamical 
systems like infinite-dimensional systems described by delay differential equa- 
tions [46], dissipative systems [44], conservative systems related to Celestial 
Mechanics problems [53|, l55] , as well as for the determination of the LCEs 
from a time series [1441 1118] . 

4 Lyapunov Characteristic Exponents: Theoretical 
treatment 

In this section we define the LCEs of various orders presenting also the basic 
theorems which guarantee their existence and provide the theoretical back- 
ground for their numerical evaluation. In our presentation we basically follow 
the fundamental papers of Oseledec [102] and of Benettin et al. [13] where all 
the theoretical results of the current section are explicitly proved. 

We consider a continuous or discrete dynamical system defined on a dif- 
ferentiable manifold S. Let ^*(x) denote the state at time t of the system 
which at time t = was at x (see equations ([3]) and ^ for the continuous 
and discrete case respectively). For the action of ^* over two successive time 
intervals t and s we have the following composition law 
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Fig. 3. Time evolution of appropriate quantities denoted by Xp ^ , P ~ 1,2,3, having 
respectively as limits for t ^ oo the first three LCEs xi, Xa, Xs, i^v two chaotic 
orbits (left panel) and one regular orbit (right panel) of the 3D Hamiltonian system 
initially studied in [34] (see [1^ for more details). In both panels Xg'^ tends to zero 
implying that Xs — 0. This is due to the fact that Hamiltonian systems have at least 
one vanishing LCE, namely the one corresponding to the direction along the flow 
(this property is explained in Section [4.5p . On the other hand, Xi a^nd X2 seem to 
get nonzero values (with xi > Xz) for chaotic orbits, while they appear to vanish 
for regular orbits (after |14J). 



The tangent space at x is mapped onto the tangent space at <?*(x) by the 
differential c?x^* according to equation ([7]). The action of <?*(x) is given by 
equation ^ for continuous systems and by equation (|12p for discrete ones. 
Thus, the action of d^'P on a particular initial deviation vector w of the 
tangent space is given by the multiplication of matrix Y(t) for continuous 
systems or Y„ for discrete systems with vector w. From equations © and 
we see that the action of d^'P* over two successive time intervals t and s 
satisfies the composition law 

d^'P'+' ^d^.^^^-p'od^-P^ (22) 

This equation can be written in the form 

R(t + s,x) = R(t,<^"(x)) • R(s,x), (23) 

where R(i, x) is the matrix corresponding to dx^^*- We note that since Y(0) = 
Yq — I2N we get dx^^w — w and R(0, x) = l2jv- A function R(i, x) satisfying 
relation (j23p is called a multiplicative cocycle with respect to the dynamical 
system 
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Let iS be a measure space with a normalized measure /x such that 

^^{S) = 1 , {<PU) = ^IiA) (24) 

for ^ C iS. Suppose also that a smooth Riemannian metric || || is defined on 
iS. We consider the multiplicative cocycle R(t, x) corresponding to dx^^* and 
we are interested in its asymptotic behavior for t —>■ ±oo. Since, as mentioned 
by Oseledec |102l . the case t +oo is analogous to the case t — > — oo, we 
restrict our treatment to the case t — )■ +oo, where time is increasing. In order 
to clarify what we are practically interested in let us consider a nonzero vector 
w of the tangent space T^S at x. Then the quantity 



At(x) 



is called the coefficient of expansion in the direction ofw. If 

lim sup - In \t (x) > 

t — ^oo t 

we say that exponential diverge occurs in the direction of w. Of course the 
basic question we have to answer is whether the characteristic exponent (also 
called characteristic exponent of order 1 ) 

lim — In A((x) 

t— >oo t 

exists. 

We will answer this question in a more general framework without re- 
stricting ourselves to multiplicative cocycles. So, the results presented in the 
following Section WA\ are valid for a general class of matrix functions, a sub- 
class of which contains the multiplicative cocycles which are of more practical 
interest to us, since they describe the time evolution of deviation vectors for 
the dynamical systems we study. 



4.1 Definitions and basic theorems 

Let At be an 71 X rt matrix function defined either on the whole real axis or 
on the set of integers, such that Aq = In, for each time t the value of function 
At is a nonsingular matrix and ||At|| the usual 2-norm of aJI. In particular, 
we consider only matrices At satisfying 

^ The 2-norm ||A|| of an n x matrix A is induced by the 2-norm of vectors, 
i. e. the usual Euclidean norm ||x|| = (X^T^i by 



and is equal to the largest eigenvalue of matrix vA^A. 
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max{||A,||,||A-i||}<e=* (25) 
with c > a suitable constant. 

Definition 4. Considering a matrix function Af as above and a nonzero vec- 
tor w of the Euclidian space the quantity 

X(^t,w) = limsupiln||Atw|| (26) 

t — ^oo t 

is called the 1-dimensional Lyapunov Characteristic Exponent or the 
Lyapunov Characteristic Exponent of order 1 (1-LCE) of At with re- 
spect to vector w. 

For simplicity we will usually refer to 1-LCEs as LCEs. 

We note that the value of the norm ||w|| does not influence the value 
of x(Af,w). For example, considering a vector /3w, with e R a nonzero 
constant, instead of w in Definition [4l we get the extra term In |/3|/t (with | | 
denoting the absolute value) in equation (|26p whose limiting value for t oo 
is zero and thus does not change the value of x(Ai, w). More importantly, the 
value of the LCE is independent of the norm appearing in equation (|26p . This 
can be easily seen as follows: Let us consider a second norm || ||' satisfying 
the inequality 

/3i||w|| < ||w||' </32||w|| 

for some positive real numbers /3i, (32 , and for all vectors w. Such norms 
are called equivalent (see e.g. |731 §5.4.7]). Then, by the above-mentioned 
argument it is easily seen that the use of norm || ||' in (^5]) leaves unchanged 
the value of x(At,w). Since all norms of finite dimensional vector spaces are 
equivalent, we conclude that the LCEs do not depend on the chosen norm. 

Let Wi, j = 1, 2, . . . ,p be a set of linearly independent vectors in R", E'p 
be the subspace generated by all and volp(At, i?^') be the volume of the 
p-parallelogram having as edges the p vectors AtWi. This volume is computed 
as the norm of the wedge product of these vectors (see Appendix for the 
definition of the wedge product and the actual evaluation of the volume) 

vo\p{At,EP) = IjAtv^fi AAtW2 A--- AAfV^fpll. 

Let also volp(Ao, E^) be the volume of the initial p-parallelogram defined by 
all Wi, since Ag is the identity matrix. Then the quantity 

^'^"^ ^ - vol^Ao, EP) 

is called the coefficient of expansion in the direction of E^ and it depends 
only on E^ and not on the choice of the linearly independent set of vectors. 
Obviously for an 1-dimensional subspace E^ the coefficient of expansion is 
||AtWi||/||wi||. If the limit 
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lim - In Xt{EP) 

t — ^oo t 

exits it is called the characteristic exponent of order p in the direction of E^. 

Definition 5. Considering the linearly independent set Wi, i ~ 1,2, . . . ,p and 
the corresponding sub space ofW^ as above, the p-dimensional Lyapunov 
Characteristic Exponent or the Lyapunov Characteristic Exponent of 
order p (p-LCE) of At with respect to subspace E^ is defined as 

x(At,^P) =limsupilnvolp(At,^P). (27) 

t — >oo t- 



Similarly to the case of the 1-LCE, the value of the initial volume volp(Ao, i?^), 
as well as the used norm, do not influence the value of x(-A-t, E^). 

From (|25l) and the Hadamard inequality (see e. g. |102j ). according to 
which the Euclidean volume of a p-parallelogram does not exceed the product 
of the lengths of its sides, we conclude that the LCEs of equations (|26|) and 
([27|) are finite. 

From the definition of the LCE it follows that 

X(At,ciWi + C2W2) < max{x(At, wi),x(At, W2)} 

for any two vectors wi,W2 G K" and ci,C2 G K with ci,C2 ^ 0, while the 
Hadamard inequality implies that if w.;, z = 1, 2, . . . , n is a basis of M" then 




where det At is the determinant of matrix Af. 

It can be shown that for any r e R the set of vectors {w G R" : x(At, Vif) < r} 
is a vector subspace of R" and that the function x( At , w) with w G R" , w 7^ 
takes at most n different values, say 

vi > V2 > ■ ■ ■ > Vs with 1 < s < n. (29) 

For the subspaces 

= {wgR" :x(At,w) <j.a (30) 

we have 

dof 

R" = D L2 D • • • D D Ls+i = {0} , (31) 

with Li^i ^ Li and x(At,w) = i>i if and only if w G \ £^+1 for i = 
1, 2, . . . , s. So in descending order each LCE 'lives' in a space of dimensionality 
less than that of the preceding exponent. Such a structure of linear spaces with 
decreasing dimension, each containing the following one, is called a filtration. 
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Definition 6. A basis Wi, i — 1,2, ... ,n o/R" is called normal ifJ27=i x(^ti w^) 
attains a minimum at this basis. In other words, the basis w^, is a normal 
basis if 

n n 
1=1 i=l 

where gi, i — 1,2, ... ,n is any other basis ofM". 

A normal basis w^, i = 1,2, ... ,n is not unique but the numbers x(At,Wi) 
depend only on At and not on the particular normal basis and are called the 
LCEs of function At- By a possible permutation of the vectors of a given 
normal basis we can always assume that x{At,Wi) > x(-A-t,W2) > ••• > 
x(At,w„). 

Definition 7. Let w^, i = 1,2, . . . ,n be a normal basis o/M" and Xi ^ X2 > 
••■ > Xn, with Xi = x(-^tjWi), i — 1,2,..., 71, the LCEs of these vectors. 
Assume that value Vi, i = l,2,...,s appears exactly ki = ki[vi) > times 
among these numbers. Then ki is called the multiplicity of value i>i and the 
collection [vi, ki) i — 1,2, . . . , s is called the spectrum of LCEs. 

In order to clarify the used notation we stress that Xij * ~ 1,2, ... ,n are 
the n (possibly nondistinct) LCEs, satisfying Xi X2 > • • • > Xm while Vi, 
i — 1,2, ... ,s represent the s {1 < s < n), different values the LCEs have, 
with vi > 1^2 > ■ ■ ■ > i^s- 

Definition 8. The matrix function At is called regular as i —> oo if for each 
normal basis Wi, i ^ 1,2, ... ,n it holds that 

n ^ 

V'xC^tiWi) = liminf-ln|det At|, 

^ ' t~tQC t 

i=l 

which, due to i28\) leads to 

lim inf — In | det j4t | = lim sup — In | det At \ 

t^oo t oo t 

guaranteeing that the limit 

hm - In I det At \ 

t^oo t 

exists, is finite and equal to 

^ n s 

lim -ln|det At| = V^xl^t^Wi) = V'fcifi. 

t—^oo t ^ — ^ ^ — ^ 

i=l 1=1 



We can now state a very important theorem for the LCEs: 
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Theorem 1. If the matrix function At is regular then the LCEs of all orders 
are given by equations i26\) and {21^ where the lim sup is substituted by lim 

t-»oo 

x(At,w) = lim iln||Atw|| (32) 

x{At,EP) = lim - \nYo\{AuEP). (33) 
In particular, for any p~dimensional subspace C M" we have 

x{At,En^J2x^r (34) 
with a suitable sequence 1 < ii < 42 < • • ■ < *p < JT-- 

The part of the theorem concerning equations ((32)l and ([33]) was proved by 
Oseledec in |102) . while equation (|34p . although was not explicitly proved in 
|102| . can be considered as a rather easily proven byproduct of the results 
presented there. Actually, the validity of equation (|34p was shown in [T3] . 

4.2 Computing LCEs of order 1 

Let us now discuss how we can use Theorem [T] for the numerical computation 
of LCEs, starting with the computation of LCEs of order 1. 

As we have already mentioned in the LCE takes at most n different 
values Vi , i — 1,2,.. .,s, 1 < s < n. If we could know a priori the sequence 
(PT|) of subspaces Li i — 1,2, s of K" we would, in principle, be able to 
compute the values of all LCEs. This could be done by taking an initial 
vector w,; E Li\ L,;+i and compute 

= lim i In j|AtWi|| , i = l,2,...,s. (35) 

t^oo t 

Now apart from Li = R" all the remaining subspaces Li, i = 2,3,. ..,s 
have positive codimension codim(Li) {— dimE" — dimL^ > 0) and thus, 
vanishing Lebesgue measure. Then a random choice of w e M" would lead to 
the computation of xi from psp . because, in principle w will belong to Li 
and not to the subspaces Li i — 2, . . . , s. Let us consider a simple example in 
order to clarify this statement. 

Suppose that Li is the usual 3-dimensional space R"^, L2 C Li is a partic- 
ular 2-dimensional plane of R'^, e. g. the plane z = 0, L3 C L2 is a particular 
1-dimensional line e. g. the x axis (Figure [4lJ a)) and the corresponding LCEs 
are Xi > X2 > X3 with multiplicities ki = k2 = ~ I. For this case we have 
dimLi = 3, dimL2 = 2, dimL3 — 1 and codim(Li) = 0, codim(L2) = 1; 
codim(L3) = 2. Concerning the measures /i of these subspaces of M.^, it is 
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Fig. 4. (a) A schematic representation of the sequence of subspaces (|31[) where 
Li identifies with R'^, L2 C Li is represented by the xy plane and the x axis is 
considered as the final subspace L3 C ^2- (b) A random choice of a vector in 
Li = will result with probability one to a vector belonging to Li and not to L2 , 
like vector wi. Vectors W2, W3 belonging respectively to L2 \ I/3 and to L3 are not 
random since their coordinates should satisfy certain conditions. In particular, the 
z coordinate of W2 should be zero, while both the 2 and y coordinate of W3 should 
vanish. The use of wi, W2, W3 in (|35|) leads to the computation of xi, X2 and X3 
respectively. 



obvious that fJ-{L2) ~ /^(^s) = 0, since the measure of a surface or of a line in 
the 3-dimensional space M'^ is zero. 

If we randomly choose a vector w G M'^ it will belong to Li and not to L^-, 
i. e. having its z coordinate different from zero and thus, equation ([35]) would 
lead to the computation of the mLCE xi- Vector Wi in Figure |4fb) represents 
such a random choice. In order to compute xi from (|35p we should choose 
vector w not randomly but in a specific way. In particular, it should belong to 
Li but not to L3, so its z coordinate should be equal to zero. Thus this vector 
should have the form w = (w\, 1112, 0) with wi,W2 £ K, W2 0, like vector W2 
in Figure HJJb) . Our choice will become even more specific if we would like to 
compute X3 because in this case w should be of the form w = 0, 0) 7^ 
with wi G M. Vector W3 of Figure |3{b) is a choice of this kind. 

From this example it becomes evident that a random choice of vector w 
in ([55)1 will lead to the computation of the largest LCE xi with probability 
one. One more comment concerning the numerical implementation of equation 
should be added here. Even if in some special examples one could happen 
to know a priori the subspaces Li i = 1, 2, . . . , s, so that one could choose 
w G Li\ Li+i with i ^ 1 then the computational errors would eventually lead 
to the numerical computation of xi- Such an example was presented in |14j . 
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4.3 Computing LCEs of order p > 1 

Let us now turn our attention to the computation of p-LCEs with p > 1. 
Equation ([M)) of Theorem [T] actually tells us that the p-LCE xi-^t,EP) can 

take at most (^^^ distinct values, i. e. as many as all the possible sums of 

p 1-LCEs out of n are. Now, as the choice of a random vector w G K", 
or in other words, of a random 1-dimensional subspace of R" produced by 
w, leads to the computation of the maximal 1-LCE, the random choice of 
a p-dimensional subspace of R", or equivalently the random choice of p 
linearly independent vectors i = 1,2, ... ,p, leads to the computation of 
the maximal p-LCE (p-mLCE) which is equal to the sum of the p largest 
1-LCEs 

p 

x{AuEP) = Y.^^- (36) 

This relation was formulated explicitly in [111 [5] and proved in |13| but was 
implicitly contained in [102] . The practical importance of equation ([55)1 was 
also clearly explained in [119^ . Benettin et al. [13] gave a more rigorous form 
to the notion of the random choice of E^, which is essential for the derivation 
of ([55]) . by introducing a condition that subspace EP should satisfy. They 
named this condition Condition R (at random). According to Condition R a 
p-dimensional space E^ C M" is chosen at random if for all j = 2, 3, . . . , s we 
have 

dim(£;P n Lj ) = max J 0, p - V fc, I (37) 



where Lj belongs to the sequence of subspaces ([51]) and ki is the multiplicity 
of the LCE V, (Definition [T]). 

In order to clarify these issues let us consider again the example presented 
in Figure m where we have three distinct values for the 1-LCEs Xi > X2 > X3 
with multiplicities fci = ^2 = fcs = 1. In this case the 2-LCE can take one of 
the three possible values Xi + X2, X2 + X3j Xi + XSi while the 3-LCE takes 
only one possible value, namely Xi + X2 + X3- 

The computation of the 2-LCE requires the choice of two linearly indepen- 
dent vectors wi, W2 and the application of equation ([55)) . The two vectors Wi, 
W2 define a 2-dimensional plane E"^ in and x(-A-f , E"^) practically measures 
the time rate of the coefficient of expansion of the surface of the parallelogram 
having as edges the vectors AfWi, AtW2. 

By choosing the two vectors wi , W2 randomly we define a random plane E'^ 
in R'^ which intersects the subspace L2 (plane xy) along a line, i. e. dim(£^^ n 
L2) — 1 and the subspace L3 (x axis) at a point, i. e. dim(i?^ni3) = (Figure 
[3Ia))- This random choice of plane E^ satisfies Condition R (|37p and thus, 
equation ([55]) leads to the computation of the 2-mLCE, namely Xi +X2- This 
result can be also understood in the following way. Plane E^ in Figure [H^a) 
can be considered to be spanned by two vectors Wi, W2 such that wi e Li 
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Fig. 5. Possible choices of the 2-dimensional space for the computation of the 
2-LCE in the example of Figure |4j where is considered as the tangent space of 
a hypothetical dynamical system. In each panel the chosen 'plane' E'^ is drawn, as 
well as one of its possible basis constituted of vectors wi, W2. (a) a random choice of 
E^ leads to a plane intersecting L2 along line e (dim(_B^ nL2) = 1) and Li at point 
A (dim(£;^ n L3) = 0). In this case equation (133) gives x(At, -B^) = xi + X2- More 
carefully made choices of E^ (which are obviously not made at random) results to 
configurations leading to the computation of xa+Xs (b) and Xi+Xs (c) from equation 
H33p . In these cases E'^ does not satisfy Condition R (|37|l since dim(_E^ n L2) = 2, 
dim(£^ n L3) = 1 in (b) and dim(£'^ n L2) = 1, dim(£^ n L3) = 1 in (c). 
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but not in its subspace L2 and W2 £ L2 but not in its subspace L3. Then the 
expansion of wi E Li\L2 is determined by the LCE xi a-nd the expansion of 
■W2 G L2\L3 by the LCE X2- These l-dimensional expansion rates result to 
an expansion rate equal to xi + X2 for the surface defined by the two vectors. 

Other more carefully designed choices of the E"^ subspace lead to the com- 
putation of the other possible values of the 2-LCE. If for example wi G L2\L^ 
and W2 G L3 (Figure [Hb)) we have E"^ = L2 with dim(i?^ n L2) = 2 and 
dim(£'2ni3) = 1. In this case the expansion of wi is determined by the LCE 
X2 and of W2 by X3, and so the computed 2-LCE is X2 + X3- Finally, a choice of 
E^ of the form presented in Figure[5Uc) leads to the computation of Xi +X3- In 
this case the plane E^ is defined by wi G ii \ ^2 and W2 G £3 and intersects 
subspaces L2 and L3 along the line corresponding to L3, i. e. dim(i?^nL2) = 1 
and dim(£^^ 0^3) = 1. It can be easily checked that for the last two choices of 
E'^ (Figures [5]Jb) and (c)), for which the computed 2-LCE does not take its 
maximal possible value, Condition R ([57)1 is not satisfied, as one should have 
expected from the fact that these choices correspond to carefully designed 
configurations and not to a random process. 

Similarly to the case of the computation of the 1-LCEs we note that, 
even if in some exceptional case one could know a priori the subspaces Li 
i = 1, 2, . . . , s, so that one could choose i = 1, 2, . . . ,p to span a particular 
subspace E^ in order to compute a specific value of the p-LCE, smaller than 
Sr=i Xi (like in Figures [Sljb) and (c)), the inevitable computational errors 
would eventually lead to the numerical computation of the maximal possible 
value of the p-LCE. 

Summarizing we point out that the practical implementation of Theorem[T] 
guarantees that a random choice of p initial vectors Wii = l,2,...,p with 1 < 
p < n generates a space E^ which satisfies Condition R ((37)) and leads to the 
actual computation of the corresponding p-mLCE, namely Xi + X2 + • ■ • + Xp- 
This statement, which was originally presented in [ITJ[3], led to the standard 
algorithm for the computation of all LCEs presented in [T3]. This algorithm 
is analyzed in Section 16.11 

4.4 The Multiplicative Ergodic Theorem 

After presenting results concerning the existence and the computation of the 
LCEs of all orders for a general matrix function Af , let us restrict our study to 
the case of multiplicative cocycles R(t, x), which are matrix functions satisfy- 
ing equation The multiplicative cocycles arise naturally in discrete and 
continuous dynamical systems as was explained in the beginning of Section |4l 
In particular, we consider the multiplicative cocycle dx^* which maps the 
tangent space at x G 5 to the tangent space at ^*(x) G S for a dynamical 
system defined on the differentiable manifold S. We recall that 5 is a measure 
space with a normalized measure /i and that is a diffeomorphism on S, 
i. e. is a measurable bijection of S which preserves the measure /i (p4)) 
and whose inverse is also measurable. We remark that in measure theory we 
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disregard sets of measure 0. In this sense is called measurable if it becomes 
measurable upon disregarding from S a set of measure 0. Quite often we will 
us the expression 'for almost all x with respect to measure fj, ' for the validity 
of a statement, implying that the statement is true for all points x with the 
possible exception of a set of points with measure 0. 

A basic property of the multiplicative cocycles is their regularity, since 
Theorem[T]guarantees the existence of characteristic exponents and the finite- 
ness of the LCEs of all orders for regular multiplicative cocycles. Thus, it is 
important to determine specific conditions that multiplicative cocycles should 
fulfill in order to be regular. Such conditions were first provided by Oseledec 
|102| who also formulated and proved the so-called Multiplicative Ergodic 
Theorem (MET), which is often referred as Oseledec's theorem. 

The MET gives information about the dynamical structure of a multiplica- 
tive cocycle R(i, x) and its asymptotic behavior for t ^ oo. The application of 
the MET for the particular multiplicative cocycle dx^?* provides the theoret- 
ical framework for the computation of the LCEs for dynamical systems. The 
MET is one of the milestones in the study of ergodic properties of dynamical 
systems and it can be considered as a sort of a spectral theorem for random 
matrix products [113j . As a testimony to the importance of this theorem one 
can find several alternative proofs for it in the literature. The original proof of 
Oseledec [1021 applies both to continuous and discrete systems. In view to the 
application to algebraic groups, Raghunathan [108j devised a simple proof of 
the MET, which nevertheless could not guarantee the finiteness of all LCEs. 
Although Raghunathan's results apply only to maps, an extension to flows, 
following the ideas of Oseledec, was given by Ruelle [114j . Benettin et al. [13] 
proved a somewhat different version of the theorem being mainly interested to 
its application on Hamiltonian flows and symplectic maps. Alternative proofs 
can also be found in [76l I141j . 

In [102j Oseledec proved that a multiplicative cocycle R(t, x) is regular 
and thus, the MET is applicable to it, if it satisfies the condition 

sup ln+ ||R±(<,x)|| e 2.1(5, Ai), 1 (38) 
t|<i 

where In"*" a = max{0. In a}. From (j38p we obtain the estimate 

||R(i,x)|| < e^("'l*l (39) 

for t ±oo for almost all x with respect to /i, where J(x) is a measurable 
function. From ([BQ]) it follows that R(t,x), considered as a function of t for 

^ We recall that a measurable function / : 5 ^ R (or C) of the measure space 
(5,/i) belongs to the space L^{S,n) if its absolute value has a finite Lebesgue 
integral, i. e. 

J \,f\dii < oo . 
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fixed X, satisfies equation P5|) . Benettin et al. [T3] considered a slightly dif- 
ferent version of the MET with respect to the one presented in [102j . Their 
version was adapted to the framework of a continuous or discrete dynamical 
system with being a diffeomorphism of class C^, i. e. both and its inverse 
are continuously differentiable. They formulated the MET for the particular 
multiplicative cocycle d^'P , which they proved to be regular. Since our presen- 
tation is mainly focused on autonomous Hamiltonian systems and symplectic 
maps we will also state the MET for the specific cocycle cix^*- The version of 
the MET we present is mainly based on [1021 11141 [T5] and combines different 
formulations of the theorem given by various authors over the years. 

Theorem 2 (Multiplicative Ergodic Theorem MET). Consider a dy- 
namical system as follows: Let its phase space S be an n-dimensional compact 
manifold with a normalized measure fi, /i(iS) — 1 and a smooth Riemannian 
metric \\ \\. Consider also a measure-preserving diffeomorphism <P* of class 
satisfying 

with t denoting time and having real (continuous system) or integer (discrete 
system) values. Then for almost all x G 5, with respect to measure fj, we have: 

1. The family of multiplicative cocycles c?x^^* : T^x.S T^pt^^^S, where T^S 
denotes the tangent space of S at point x, is regular. 

2. The LCEs of all orders exist and are independent of the choice of the 
Riemannian metric of S. 

In particular, for any w G T^S the finite limit 

X(x,vir)= lim -ln\\d^<P*-w\\ (40) 

t^oo t 

exists and defines the LCE of order 1 (1-LCE). There exists at least one 
normal basis Vj, i — 1, 2, . . . , rt of T^S for which the corresponding (pos- 
sibly nondistinct) 1-LCEs Xi(x) — x(x, Vj) are ordered as 

Xi(x) >X2(x) > ••• >x«(x). (41) 

Assume that the value ^'i(x), i = 1, 2, . . . , s with s — s(x), 1 < s < n ap- 
pears exactly fci(x) — ki^K-^Vi) > times among these numbers. Then the 
spectrum of LCEs (z/i(x), ki^x.)), i — 1,2,. ..,s is a measurable function 
of X, and as w ^ varies in T^S, x(x,w) takes one of these s different 
values 

1^1 (x) > j^2(x) > • • • > J^s(x). (42) 

R also holds 

^ 1 
Vfcj(x)zyj(x) = lim -ln|detdx'?>*|. (43) 

■'- — ' t^oo t 

i=l 

For any p-dimensional (1 < p < n) subspace E^ C 7^5, generated by a 
linearly independent set Wi, i — 1,2, . . . ,p the finite limit 
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x{x,EP)^ lim -lnvoUd^<P\EP), (44) 

t^OG t 

where volp(dx^*, -E^) is the volume of the p-parallelogram having as edges 
the vectors d^^^'Wi, exists and defines the LCE of order p (p-LCE). The 
value of x{^, E^) is equal to the sum of p 1-LCEs Xi(x), i = 1, 2, . . . , n. 
3. The set of vectors 

L,(x) = {w e T^S : x(x,w) < v,{-k)} , l<i<s 

is a linear subspace of T^S satisfying 

def 

T^S= ii(x) D L2(x) D • • • D L,(x) D i,+i(x) J {0}. (45) 

//w £ Li{x) \ii+i(x) then x(x, w) = t'i(x) for i = 1,2, ... ,s. The multi- 
plicity fci(x) of values ^i(x) is given by fci(x) = diniLi(x) — dimii+i(x). 
-i^. r/ie symmetric positive-defined matrix 

/ „ N l/2t 

A^^ lim (rT(i). F(0) 

exists. Y(t) is the matrix corresponding to d^^^ and is defined by equa- 
tions nop and I113\) for continuous and discrete dynamical systems re- 
spectively. The logarithms of the eigenvalues of /Ix are the s distinct 1- 
LCEs \4-2^ of the dynamical system. The corresponding eigenvectors are 
orthogonal (since A^x. is symmetric) , and for the corresponding eigenspaces 
Vi(x), V2(x), . . . , "I4(x) we have 

s 

fc,(x) = diml/,(x) , L,(x)=01/,(x) /or z = l,2,...,s. 

Thus, Tx_S is decomposed as 

T^S = 14 (x) © l^2(x) © • • ■ ® V;(x), 
and for every nonzero vector w £ V^i(x), i — 1,2, . . . ,s, we get 

X(x,w) = i/,(x). 

A short remark is necessary here. The regularity of dx^*, which guarantees 
the validity of equations (|40|) and (|44)) and the finiteness of the LCEs of 
all orders, should not be confused with the regular nature of orbits of the 
dynamical system. Regular orbits have all their LCEs equal to zero (see also 
the discussion in Section [O]) . 

Benettin et al. [TTl [T3] have formulated also the following theorem which 
provides the theoretical background for the numerical algorithm they pre- 
sented in [T^ for the computation of all LCEs. 
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Theorem 3. Under the assumptions of the MET, the p-LCE of any p- 
dimensional subspace E^ C Ty^S satisfying Condition R ^37^ , is equal to the 
sum of the p largest l~LCEs 

xi^.EP) lim -\nvolp{dy<P\EP) = Vx.W- (46) 

t^OG t ^ ^ 

i=l 



4.5 Properties of the spectrum of LCEs 

Let us now turn our attention to the structure of the spectrum of LCEs for 
A^D autonomous Hamiltonian systems and for 2Nd symplectic maps, which 
are the main dynamical systems we are interested in. Such systems preserve 
the phase-space volume, and thus, the r. h. s. of ([151) vanishes. So for the sum 
of all the 1-LCEs we have 

2N 

J2M^) = 0- (47) 

The symplectic nature of these systems gives indeed more. It has been 
proved in [T3] that the spectrum of LCEs consists of pairs of values having 
opposite signs 

X^{x) = ~X2N-^+l{^) , l=l,2,...,N. (48) 

Thus, the spectrum of LCEs becomes 

Xl(x) > X2(x) > • • • > Xw(x) > -XAr(x) > • • • > -X2(x) > -Xl(x). 

For autonomous Hamiltonian flows we can say something more. Let us 
first recall that for a general differentiable flow on a compact manifold without 
stationary points at least one LCE must vanish [T31 [7D] . This follows from the 
fact that, in the direction along the flow a deviation vector grows only linearly 
in time. So, in the case of a Hamiltonian flow, due to the symmetry of the 
spectrum of LCEs ([48]), at least two LCEs vanish, i. e. 

XJv(x) = X7v+i(x) = 0, 

while the presence of any additional independent integral of motion leads to 
the vanishing of another pair of LCEs. 

Let us now study the particular case of a periodic orbit of period T, such 
that ^-^(x) = X, following [9l[T2]. In this case c?x^"^ is a linear operator on 
the tangent space 7^5 so that for any deviation vector w{Q) G 7^5 we have 

w(r) = Y-w(0), (49) 

where Y is the constant matrix corresponding to d^ip^ . Suppose that Y has 
2N (possibly complex) eigenvalues A^, i — 1,2, ... , 2N, whose magnitudes can 
be ordered as 
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|Ai| > IA2I > ... > \X2n\- 

Let Wi, i — 1,2, ... , 2N, denote the corresponding unitary eigenvectors. Then 
for w(0) — "Wi equation (j49p imphes 

w(fcr) = Aj^w, , fc = 1,2,... (50) 

and so we conclude from pO)) that 

X(x,w,) = L\n\X,\=x^{^), i = l,2,...,2N. 

Furthermore for a deviation vector 

w(0) = CiWi + C2W2 + . . . + C2N^2N 

with Ci & i — 1,2, ... , 2N , it foUows from (|50p that the first nonvanishing 
coefficient Ci eventuaUy dominates the evolution of w(t) and we get x(x, w) = 
Xi. In this case we can define a filtration similar to the one presented in 
dUl) by defining Li = [wi, W2, . . . , W2Ar] = T^S, L2 = [w2, . . . , W2Ar], 
L2N — [w27v], L2N+1 = [0], where [ ] denotes the linear space spanned by 
vectors Wi,W2, . . . , W27v and so on. It becomes evident that a random choice 
of an initial deviation vector w(0) G T^S will lead to the computation of the 
mLCE Xi(x) since, in general, w(0) & Li\ L2. 

So, in the case of an unstable periodic orbit where |Ai| > 1 we get 
Xi(x) > 0, which implies that nearby orbits diverge exponentially from the 
periodic one. These orbits are not called chaotic, although their mLCE is 
larger than zero, but simply 'unstable'. In fact, unstable periodic orbits exist 
also in integrable systems. Since the measure of periodic orbits in a general 
dynamical system has zero measure, periodic orbits (stable and unstable) are 
rather exceptional. 

In the general case of a nonperiodic orbit we are no more allowed to use 
concepts as eigenvectors and eigenvalues because the linear operator dx^* 
maps T^iS into '^(p^f^^^'S 7^ %cS, while eigenvectors are intrinsically defined 
only for linear operators of a linear space into itself. Nevertheless, in the case of 
nonperiodic orbits the MET proves the existence of the LCEs and of filtration 
(|45|) . In a way, the MET provides an extention of the linear stability analysis 
of periodic orbits to the case of nonperiodic ones, although one should always 
keep in mind that the LCEs are related to the real and positive eigenvalues 
of the symmetric, positive-defined matrix Y'^{t) ■ Y{t) [S31IHH]. On the other 
hand, linear stability analysis involves the computation of the eigenvalues 
of the nonsymmetric matrix Y(i), which solves the linearized equations of 
motion (jlOp for Hamiltonian flows or (fT3|) for maps. These eigenvalues are 
real or come in pairs of complex conjugate pairs and, in general, they are not 
directly related to the LCEs which are real numbers. 

An important property of the LCEs is that they are constant in a connected 
chaotic domain. This is due to the fact that every nonperiodic orbit in the 
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same connected chaotic domain covers densely this domain, thus, two different 
orbits of the same domain are in a sense dynamicaUy equivalent. The unstable 
periodic orbits in this chaotic domain have in general LCEs that are different 
from the constant LCEs of the nonperiodic orbits. This is due to the fact 
that the periodic orbits do not visit the whole domain, thus, they cannot 
characterize its dynamical behavior. In fact, different periodic orbits have 
different LCEs. 



5 The maximal LCE 

From this point on, in order to simplify our notation, we will not explicitly 
write the dependence of the LCEs on the specific point x S 5. So, in practice, 
considering that we are referring to a specific point x G 5, we denote by Xi 

(v) 

the LCEs of order 1 and by Xi the LCEs of order p. 

For the practical determination of the chaotic nature of orbits a numerical 
computation of the mLCE xi can be employed. If the studied orbit is regular 
Xi = 0, while if it is chaotic xi > Oi implying exponential divergence of nearby 
orbits. The computation of the mLCE has been used extensively as a chaos 
indicator after the introduction of numerical algorithms for the determination 
of its value at late 70's [H [Ml li 131 H] • 

Apart from using the mLCE as a criterion for the chaoticity or the regular- 
ity of an orbit its value also attains a 'physical' meaning and defines a specific 
time scale for the considered dynamical system. In particular, the inverse of 
the mLCE, which is called Lyapunov time 

tL^— (51) 

gives an estimate of the time needed for a dynamical system to become chaotic 
and in practice measures the time needed for nearby orbits of the system to 
diverge by e (see e. g [30l p. 508]). 

5.1 Computation of the mLCE 

The mLCE can be computed by the numerical implementation of equation 
(PP)) . In Section 14.21 we showed that a random choice of the initial deviation 
vector w(0) G %t.S leads to the numerical computation of the mLCE. We 
recall that the deviation vector w(t) at time t > is determined by the action 
of the operator c?x^ on the initial deviation vector w(0) according to equation 

m 

w{t) = d^$*vf{0). (52) 

This equation represents the solution of the variational equations ([5]) or the 
evolution of a deviation vector under the action of the tangent map pT|) . 
and takes the form ^ and IT^ respectively. We emphasize that, both the 
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variational equations and the equations of the tangent map are hnear with 
respect to the tangent vector w, i. e. 

dx??* (aw) = adx?^*w, for any a £ M. (53) 

In order to evaluate the mLCE of an orbit with initial condition x(0), one 
has to follow simultaneously the time evolution of the orbit itself and of a 
deviation vector w from this orbit with initial condition w(0). In the case of 
a Hamiltonian flow (continuous time) we solve simultaneously the Hamilton 
equations of motion ^ for the time evolution of the orbit and the variational 
equations ([5]) for the time evolution of the deviation vector. In the case of a 
symplectic map (discrete time) we iterate the map (jU for the evolution of 
the orbit simultaneously with the tangent map (jlip . which determines the 
evolution of the tangent vector. The mLCE is then computed as the limit for 
f — *■ oo of the quantity 

, \Kio)^'-^m 1 , iiw(t)ii 

^'^'^ 1 ||w(o)|| " 7 Mol' ^ ' 

often called finite time mLCE. So, we have 

XI = lim Xi{t). (55) 

t — *oo 

The direct numerical implementation of equations ([54]) and ([55]) for the 
evaluation of xi meets a severe difficulty. If, for example, the orbit under study 
is chaotic, the norm l|w(t)|| increases exponentially with increasing time t, 
leading to numerical overflow, i. e. ||w(t)|| attains very fast extremely large 
values that cannot be represented in the computer. This difficulty can be 
overcome by a procedure which takes advantage of the linearity of dx^* (|53|) 
and of the composition law (j22p . Fixing a small time interval r we express 
time t with respect to r as t = kr, k = 1,2,.... Then for the quantity Xi{t) 
we have 

^^(,,) 1 l,.l|w(fcr)ll 



kr ||w(0)|| 

1^^^^ ||w(fcr)|| ||w((fc-l)r)|| ||w(2r)|| ||w(r)|| 



kr Vl|w((fc-l)r)||||w((fc-2)r)|| ||w(r)|| ||w(0)| 



1 



\v^(iT) 



krj^, ||w((i-l)r)|| 



i=l 
k 



Denoting by Z?o the norm of the initial deviation vector w(0) 



^0- l|w(0)||. 
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we get for the evolved deviation vector at time t = kr 



dx(0)<?'" W(0) '^<?<'-^'-(x(0))'^ w(0)|| 



||dx(o)<?**"''"w(0)|| ^0 
Let us now denote by 

U ^ rfx(0)g'^'-^^^w(0) ^ 

\Kio)<i>^^-'^^M0)\\ 

the deviation vector at point <l'^^~^^'^ (x(0)) having the same direction with 
w((z — l)r) and norm Do, and by Di its norm after its evolution for r time 
units 

A = |IV-)-(x(o))'^'^((*-l)^)ll- 
Using this notation we derive from equation (j57[) 



1 Nx(o)<^''^w(0)|| A , 

In = ln— =lnQ;j, (58) 

Nx(o)'^^^"''^w(0)|| Do 

with ai being the local coefficient of expansion of the deviation vector for a 
time interval of length r when the corresponding orbit evolves from position 
(^(»-i)'^(x(0)) to position ^''^(x(O)) (Inai/r is also called stretching number 
[135] [301 p. 257]). 

From equations ([55]) . and ([SB)) we conclude that the mLCE xi can be 
computed as 

1 ^ D. 1 ^ 

Xi = lim Xi{kT) = lim - — > In — - = lim - — > Inoii. (59) 

fe^oo k^oc KT ^ — ' Do k~*oo KT ^ — ' 

i—1 i—1 

Since the initial norm Dq can have any arbitrary value, one usually set it to 
Do = 1. Equation (|59p implies that practically xi is the limit value, for t oo, 
of the mean of the stretching numbers along the studied orbit [HI [571 1135j . 



5.2 The numerical algorithm 

In practice, for the evaluation of the mLCE we follow the evolution of a unitary 
initial deviation vector w(0) = w(0), l|w(0)|| = Do = 1 and every t — t time 
units we replace the evolved vector w(fcr), k = 1,2,..., by vector w(fcT) 
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having the same direction but norm equal to 1 (||w(fcr)|| = 1). Before each 
new renormalization the corresponding is computed and Xi is estimated 
from equation (j59p . 

More precisely at t = r we have ai — ||w(t)||. Then we define a unitary 
vector w(r) by renormalizing w(t) and using it as an initial deviation vector 
we evolve it along the orbit from x(t) to x(2t) according to equation ([52|) . 
having w(2r) = dx(T)^^ "^(t). Then we define a2 — ||w(2t)|| and we estimate 
Xi (see Figure [6]) . We iteratively apply the above described procedure until 



Fig. 6. Numerical scheme for the computation of the mLCE xi- The unitary 
deviation vector w((i — l)r), i = 1,2,..., is evolved according to the variational 
equations (|8]) (continuous time) or the equations of the tangent map (discrete 
time) for t — t time units. The evolved vector w(ir) is replaced by a unitary 
vector w(jr) having the same direction with w(jr). For each successive time interval 
[{i — l)r, ir] the quantity a, — ||w(jr)|| is computed and xi is estimated from 
equation lf59)) . 

a good approximation of xi is achieved. The algorithm for the evaluation of 
the mLCE xi is described in pseudo-code in Table [1] 

Instead of utilizing the variational equations or the tangent map for the 
evolution of a deviation vector in the above described algorithm, one could 
integrate equations ([2|) or iterate equations ([4]) for two orbits starting nearby 
and estimate w(i) by difference. Indeed, this approach, infiuenced by the rough 
idea of divergence of nearby orbits introduced in [75] , was initially adopted for 
the computation of the mLCE [101 IMl IH] ■ This technique was abandoned after 
a while as it was realized that the use of explicit equations for the evolution 



w(2t) 



w(0)=w(0) 




x(0) 
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Table 1. The algorithm for the computation of the mLCE xi as the limit for t ^ oo 
of Xi{t) according to equation H59|) . The program computes the evolution of Xi{t) 
as a function of time f up to a given upper value of time t = Tm or until X\{t) 
attains a very small value, smaller than a low threshold value Xim- 

Input: 1. Hamilton equations of motion ([2]) and variational equations ([HI, or 
equations of the map Q and of the tangent map 

2. Initial condition for the orbit x(0). 

3. Initial unitary deviation vector w(0). 

4. Renormalization time r. 

5. Maximal time: Tm and minimum allowed value of Xi{t): X\ra- 
Step 1 Set the stopping flag, SF ^ 0, and the counter, <— 1. 

Step 2 While (SF = 0) Do 

Evolve the orbit and the deviation vector from time t = (fc — l)r 

to t = fcr, i. e. Compute x(fcr) and w(fcT). 
Step 3 Compute current value of ~ ||w(A;r)||. 

Compute and Store current value of Xi{kT) = -j^ X^iLi ludi- 
Step 4 Renormalize deviation vector by Setting w(fcT) ^ ■wikr) / au- 
Step 5 Set the counter fc <— + 1. 

Step 6 If [(fcr > Tm) or (Xi((fc - l)r) < Then 
Set SF^ 1. 
End If 
End While 
Step 7 Report the time evolution of X\{t). 



of deviation vectors was more reliable and efficient [34l I119[ [T4] , although in 
some cases it is used also nowadays (see e. g. [145] ). 

5.3 Behavior of for regular and chaotic orbits 

Let us now discuss in more detail the behavior of the computational scheme 
for the evaluation of the mLCE for the cases of regular and chaotic orbits. 

The LCE of regular orbits vanish [IQj .23J due to the linear increase with 
time of the norm of deviation vectors. We illustrate this behavior in the case 
of an A^D Hamiltonian system, but a similar analysis can be easily carried out 
for symplectic maps. In such systems regular orbits lie on A^-dimensional tori. 
If such tori are found around a stable periodic orbit, they can be accurately 
described by N formal integrals of motion in involution, so that the system 
would appear locally integrable. This means that we could perform a local 
transformation to action-angle variables, considering as actions Ji, J2, . . . , Jn 
the values of the N formal integrals, so that Hamilton's equations of motion, 
locally attain the form 

j, = 0, =cc;,(Ji, J2,..., Jat), z = l,2,...,iV. (60) 

These equations can be easily integrated to give 



Ji{t) — Jio, 6i{t) — 9if) + uji{Ji(), J20, . . . , Jno) t, i — 1,2, . . . , N, 
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where J^o: ^iO: * — 1, 2, . . . , iV are the initial conditions of the studied orbit. 

By denoting as ^i, Ty^, z = 1, 2, . . . , iV small deviations of Ji and Bi respec- 
tively, the variational equations ([5]) of system ([SO)) , describing the evolution 
of a deviation vector are 



JV 



ii = 0, i]i = ^ • 



1,2,. 



where 



9wj 

a7" 



, = 1,2,...,7V, 



and Jq = (Jioj >/20i ■ • ■ i "^Aro) = constant, represents the A'^-dimensional vector 
of the initial actions. The solution of these equations is: 



77.W = ^.(o)+[Ef=i^y-0(o) 



1,2,. 



(61) 



From equations (pT|) we see that an initial deviation vector w(0) with coor- 
dinates fi(0), z = 1, 2, . . . , TV in the action variables and r7i(0), i = 1, 2, . . . , iV 
in the angles, i. e. w(0) = (6(0), 6(0), ■ • ■ , CjvW, mW, ?72(0), . . . , w(0)), 
evolves in time in such a way that its action coordinates remain constant, 
while its angle coordinates increase linearly in time. This behavior implies an 
almost linear increase of the norm of the deviation vector. To see this, let us 
assume that vector w(0) has initially unit magnitude, i. e. 



N 



N 



whence the time evolution of its norm is given by 

2' 



l|w(i)l|-U + 



N / N 

i=i \j=i 



N 



N 



2E r^'WE'^^.c^w 



1/2 



This implies that the norm for long times grows linearly with t 

||w(0|| ext. 



(62) 



So, from equation ((54l) we see that for long times Xi{t) is of the order 0{\nt/t), 
which means that Xi{t) tends asymptotically to zero, as t — > oo like t~^. This 
asymptotic behavior is evident in numerical computations of the mLCE of 
regular orbits, as we can see for example in the left panel of Figured) 

The asymptotic behavior of Xi{t) for regular orbits, described above, rep- 
resents a particular case of a more general estimation presented in [63]. In 
particular, Goldhirsch et al. [S3] showed that, in general, after some initial 
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transient time the value of the niLCE xi is related to its finite time estima- 
tion by 

X,W=Xi + ^, (63) 

where 6 is a constant and z{t) is a 'noise' term of zero mean. According to their 
analysis, this approximate formula is valid both for regular and chaotic orbits. 
It is easily seen that from (|63|) we retrieve again the asymptotic behavior 
Xi{t) oc t^^ for the case of regular orbits (xi = 0). 

In the case of chaotic orbits the variation of Xi (t) is usually irregular for 
relatively small t and only for large t the value of Xi (t) stabilizes and tends 
to a constant positive value which is the mLCE xi- If for example the value of 
Xi is very small then initially, for small and intermediate values of t, the term 
proportional to t"-^ dominates the r. h. s. of equation and Xi{t) oc t~^. 
As t grows the significance of term (b + z{t))/t diminishes and eventually the 
value of xi becomes dominant and Xi{t) stabilizes. It becomes evident that 
for smaller values of xi the larger is the time required for Xi{t) to reach 
its limiting value, and consequently Xi{t) behaves as in the case of regular 
orbits, i. e. Xi{t) oc t^^ for larger time intervals. This behavior is clearly seen 
in Figure [7] where the evolution of Xi{t) of chaotic orbits with small mLCE is 
shown. In particular, the values of the mLCE are xi ~ 8- 10^'^ (left panel) and 
Xi ~ 1-6 • 10~^ (right panel). In both panels the evolution of Xi{t) of regular 
orbits (following the power law (x t~^) is also plotted in order to facilitate the 
comparison between the two cases. 

6 Computation of the spectrum of LCEs 

While the knowledge of the mLCE xi can be used for determining the regular 
{Xi = 0) or chaotic (xi > 0) nature of orbits, the knowledge of part, or of the 
whole spectrum of LCEs, provides additional information on the underlying 
dynamics and on the statistical properties of the system, and can be used for 
measuring the fractal dimension of strange attractors in dissipative systems. 

In Section 1431 it was stated that, for Hamiltonian systems the existence of 
an integral of motion results to a pair of zero values in the spectrum of LCEs. 
As an example of such case we refer to the Hamiltonian system studied in 
[12] . This system has one more integral of motion apart from the Hamiltonian 
function and so 4 LCEs were always found to be equal to zero. Thus, the 
determination of the number of LCEs that vanish can be used as an indicator 
of the number of the independent integrals of motion that a dynamical system 
has. 

It has been also stated in Section 14.51 that the spectrum of the LCEs of 
orbits in a connected chaotic region is independent of their initial conditions. 
So, we have a strong indication that two chaotic orbits belong to connected 
chaotic regions if they exhibit the same spectrum. As an example of this 
situation we refer to the case studied in [3] of two chaotic orbits of a 16D 
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Fig. 7. Evolution of Xi{t) (denoted as Ln) with respect to the discrete time t 
(denoted as N) in log-log scale for regular (grey curves) and chaotic (black curves) 
orbits of the 4d map (|16|) (left panel) and of a 4d map composed of two coupled 2d 
standard maps (right panel) (see [122| for more details). For regular orbits Xi{t) 
tends to zero following a power law decay, Xi{t) oc t~^. For chaotic orbits Xi{t) 
exhibits for some initial time interval the same power law decay before stabilizing 
to the positive value of the mLCE xi- The length of this time interval is larger 
for smaller values of xi- The chaotic orbits have Xi ~ 8 ' 10~^ (left panel) and 
Xi ^ 1.6 • 10"'' (right panel) (after [122] ). 

Hamiltonian system having similar spectra of LCEs but very different initial 
conditions. 

Vice versa, the existence of different LCEs spectra of chaotic orbits pro- 
vides strong evidence that these orbits belong to different chaotic regions of 
the phase space that do not communicate. In |14j two chaotic orbits, pre- 
viously studied in [34j . were found to have significantly different spectra of 
LCEs and they were considered to belong to different chaotic regions which 
were called the 'big' (corresponding to the largest xi) and the 'small' chaotic 
sea. It is worth mentioning that the numerical results of [Tl] suggested the 
possible existence of an additional integral of motion for the 'small' chaotic 
sea, since X2 seemed to vanish. This assumption was in accordance to the 
results of [34| where such an integral was formally constructed. 
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The spectrum of LCEs is also related to two important quantities namely, 
the metric entropy, also called Kolmogorov-Sinai (KS) entropy h, and the 
information dimension Di, which are trying to quantify the statistical prop- 
erties of dynamical systems. For the explicit definition of these quantities, 
as well as detailed discussion of their relation to the LCEs the reader is re- 
ferred for example to [3 |46l El E] [92l p. 304-305] for the KS entropy and 
to [791 HHl Hll [Ml HI for the information dimension. 

In particular, Pesin [106j showed that under suitable smoothness condi- 
tions the relation between the KS entropy h and the LCEs is given by 



M 



Xi(x)>0 



where the sum is extended over all positive LCEs and the integral is defined 
over a specified region M. of the phase space S. 

Kaplan and Yorke [75] introduced a quantity, which they called the Lya- 
punov dimension 

Dl = J + Si^^ (64) 

where j is the largest integer for which Xi + X2 + • • ■ + Xj ^ 0- The Kaplan- 
Yorke conjecture states that the information dimension Di is equal to the 
Lyapunov dimension D^, i. e. 

Di^Dl, (65) 

for a typical system, and thus, it can be used for the determination of the 
fractal dimension of strange attractors. The meaning of the word 'typical' 
is that it is not hard to construct examples where equation (l65|) is violated 
(see e. g [H]). But the claim is that these examples are pathological in that 
the slightest arbitrary change of the system restores the applicability of (|65|) 
and that such violation has 'zero probability' of occurring in practice. The 
validity of the Kaplan- Yorke conjecture has been proved in some cases [1461 
[57] although a general proof has not been achieved yet. We note that in the 
case of a 2A^D conservative system Dl \s equal to the dimension of the whole 
space, \. e. Dl — 2N, because j = 2N in since Yl'i=i Xi — ^ according to 



equation p7)) . 

So, it becomes evident that developing an efficient algorithm for the nu- 
merical evaluation of few or of all LCEs is of great importance for the study of 
dynamical systems. In this section we present the different methods developed 
over the years for the computation of the spectrum of LCEs, focusing on the 
method suggested by Benettin et al. [14], the so-called standard method. 

6.1 The standard method for computing LCEs 

The basis for the computation of few or even of all LCEs is Theorem |3l 
which states that the computation of a p-LCE from equation (|44l) , considering 
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a random choice of p (1 < p < 2A^) linearly independent initial deviation 
vectors, leads to the evaluation of the p-mLCE Xi ; which is equal to the 



sum of the p largest 1-LCEs (|46| . 

In order to evaluate the p-mLCE of an orbit with initial condition 
x(0), one has to follow simultaneously the time evolution of the orbit it- 
self and of p linearly independent deviation vectors with initial conditions 
Wi(0), W2(0), . . . , Wp(0) (using the variational equations ([5]) or the equations 
of the tangent map (HU)). Then, the p-mLCE is computed as the limit for 
t — > cxD of the quantity 

X<-PHt) - ^ In ^"^^ (^^(o)*^* wi(0), dx(o)<?'* W2(0), ■ ■ ■ , rfx(o)'^* Wp(0)) 

^' t VOlp(wi(0),W2(0),...,Wp(0)) 

t ||wi(0)Aw2(0)A...Awp(0)|l t ||ALiW.(0)l|' ^ ' 
which is also called the finite time p-mLCE. So we have 

Xf^ ^Xi+X2 + --- + Xp = lim X^P\t). (67) 

We recall that the quantity volp (wi, W2, . . . , vifp) appearing in the above 
definition is the volume of the p-parallelogram having as edges the vectors 
wi, W2, • ■ • , vifp (see equations ()106p and (|105|) in Appendix E|. 

The direct numerical implementation of equations (|66p and (|67p faces one 
additional difficulty apart from the fast growth of the norm of deviation vec- 
tors discussed in Section 15.11 This difficulty is due to the fact that when at 
least two vectors are involved (e. g. for the computation of x\ )i the angles 
between their directions become too small for numerical computations. 

This difficulty can be overcome on the basis of the following simple remark: 
an invertible linear map, as c?x(o)^ i maps a linear p-dimensional subspace 
onto a linear subspace of the same dimension, and the coefficient of expansion 
of any p-dimensional volume under the action of any such linear map (like 
for example HAf^i / IIAr=i w,:(0)|| in our case) does not depend on the 

initial volume JJ:J. Since the numerical value of HAiLi '^i(0)ll does not depend 
on the choice of the orthonormal basis of the space (see Appendix [Al for more 
details), in order to show the validity of this remark we will consider an 
appropriate basis which will facilitate our calculations. 

In particular, let us consider an orthonormal basis {ei, §2, . . . , e^} of the 
p-dimensional space C 7^(o)5 spanned by {wi(0), W2(0), . . . ,Wp(0)}. This 
basis can be extended to an orthonormal basis of the whole 2A'^-dimensional 
space {§1, §2, . . . , ep, Gp+i, . . . , e2N} and E''' C T^fQ^S can be written as the 
direct sum of E^ and of the (2iV — p)-dimensional subspace E' spanned by 
{ep+i, . . . , e2Af} 

%,^o)S^EP®E'. 

Consider also the 2A^ x p matrix W(0) having as columns the coordinates of 
vectors Wi(0), i = 1,2, ... ,p with respect to the complete orthonormal basis 
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1,2,..., 27V, in analogy to equation (I102p . Since Wi(0) G i?^ this 



matrix has the form 



W(0) 







W(0) 

{2N-p)xp 



where W(0) is a square px p matrix and 0^2N-p)xp is the {2N —p) x p matrix 
with all its elements equal to zero. Then, according to equations (|105|) and 
(|106p the volume of the initial p-parallelogram is 



detW(O) 



(68) 



since detW (0) = detW(O) for the square matrix W(0). 

Each deviation vector is evolved according to equation ((T]) and it can be 
computed through equation ([9]) or ([12]), with Y(t) being the 2N x 2N matrix 
representing the action of (ix(o)^*- By doing a similar choice for the basis of 
the ^^*(x(o))'^ space, equation (|102p gives for the evolved vectors 

[wi(t) W2(t) •■• wp(t)] = [ei e2 ■••ep]-Y(i)-W(0)= [ei ea • ■ • e^] .W(i). 
Writing Y{t) as 

Y(t)= [Yi(t) Y,it)] 

where Yi(t) is the 2N x p matrix formed from the first p columns of Y(t) and 
Y2(t) is the 2N x {2N — p) matrix formed from the last 2N — p columns of 
Y(i), W(i) assumes the form 

W(i) = Yi(i) • W(0). 
Then from equation (|105[) we get 



W (0)-Y^(t).Yi(t).W(0) 



detW (O)dct (Y?'(i) ■ Yi(t)) detW(O) 



det W(0) I Jdet (y^ (t) • Yi (t)) . 



(69) 



Thus, from equations (|68p and ([69|l we conclude that the coefficient of 
expansion 

llALiw.(i)|| 



IIALiw.(o)| 



= Wdet 



(Y^(i).Yi(O 



does not depend on the initial volume but it is an intrinsic quantity of the 
subspaces defined by the properties of dx(o)^*- Note that in the particular 
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case of p = 27V the coefficient of expansion is equal to |detY(i)| in accor- 
dance to equation (j43p . An alternative way of expressing this property is 
that, for two sets of linearly independent vectors {vifi(O), vir2(0), . . . ,Wp{0)} 
and {f i(0), £2(0), . . . , fp(0)} spanning the same p-dimensional subspace of 
7^(0)5, the relation 



IIAr=iW.(o)|| 



(70) 



holds [nn]. 

Let us now describe the method for the actual computation of the p- 
mLCE. Similarly to the computation of the mLCE we fix a small time interval 
r and define quantity X'-P^t) ([66]) as 



k 



A^=irfx(o)<?'"w,(0)| 



X(p)(fcT) — > 'In 



1 ^ 

i—1 



(p) 



(71) 



(v) 

where 7, , i = 1,2,..., is the coefficient of expansion of a p-dimensional vol- 
ume from t — {i — l)r to t = it. According to equation ([70]) ^[^"^ can be com- 
puted as the coefficient of expansion of the p-parallelogram defined by any p 
vectors spanning the same p-dimensional space. A suitable choice for this set is 
to consider an orthonormal set of vectors {wi((i — l)r), W2((i — 1)t), . . . , Wp((z — l)r)} 
giving to equation (|7ip the simplified form 



k 



(p) 



kr ^ 

i=l 



In 



A d^{{^-l)r)^' W,((i- l)r) 



Thus, from equations (|67|) and ((72)) we get 



(72) 



(p) 
Xi 



Xi + X2 H h Xp 



1 

lim - — In 7}^^ 
fe^oo kr ^ 
1=1 



(73) 



for the computation of the p-mLCE. This equation is valid for 1 < p < 2A^ 
since in the extreme case of p = 1 it is simply reduced to equation (|59p with 
= 7,p''. In order to estimate the values of Xi, i = 1, 2, . . . ,p, which is our 
actual goal, we compute from ([75]) all the x^^'' quantities and evaluate the 
LCEs from 

xf^-xr'\ ^ = 2,3,...,p (74) 



X» 



with xl'^ = xi [ns]. 

Benettin et al. [M] noted that the p largest l~LCEs can be evaluated at 
once by computing the evolution of just p deviation vectors for a particular 
choice of the orthonormalization procedure, namely performing the Gram- 
Schmidt orthonormalization method. 
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Let us discuss the Gram-Schmidt orthonormahzation method in some de- 
tail. Let ^^{iT), i — 1, 2, ... ,p be the evolved deviation vectors Wj((i — l)r) 
from time t = {i — 1)t to t = it. From this set of linearly independent vectors 
we construct a new set of orthonormal vectors Wj (ir) from equations 



Ui(iT) = wi(ir), 7ii = ||ui(ir)||, Wi(zr) 



Ui(»t) 

7h 



U2(iT) = W2(ir) - (w2(ir), Wi(ir)}wi(ir), 72i = ||u2(ir)||, W2(ir) 
I - (w3(ir), wi(r 
73i = l|u3(iT)||, W3(iT) 



(75) 

U2(^t) 

72i 



U3(iT) = W3(ir) - (w3(ir), wi(ir))wi(ir) - (wr3(ir), W2(ir))w2(ir), 

U3(iT) 



73j 



which are repeated up to the computation of Wp(iT). We remark that (w, u) 
denotes the usual inner product of vectors vif, u. The general form of the 
above equations, which is the core of the Gram-Schmidt orthonormahzation 
method, is 



fc-i 



Ufc(ir) = viffe(iT) - ^(wrfe(ir), Wj(iT))wj(ir), 

Ufc(ir) 



Iki = ||Ufe(«T)||, Wfe(iT) = 



Iki 



(76) 



for 1 <k <p. 

As we will show in Section 16.31 the volume of the p-parallelogram having 
as edges the vectors dx((i-i)T)^^ w-j ((*"!)''') = ■Wj(ir), j = 1, 2, . . . ,p is equal 
to the volume of the p-parallelogram having as edges the vectors Uj(iT), i. e. 



/\ rfx((»-l)r)<^^Wj((i - l)r) 



(77) 



Since vectors Uj(iT) are normal to each other, the volume of their p- 
parallelogram is equal to the product of their norms. This leads to 



„(p) 



(78) 



Then, equation ([75)1 takes the form 
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Using now equation (ffi)) we are able to evaluate the 1-LCE Xp as 

(p) (p-i) 1 rij^i 1 v-^ 

In conclusion we see that the value of the 1-LCE Xp with 1 < p < 2N can 
be computed as the limiting value, for t —^ oo, of the quantity 

1 ^ 

i=l 

1. 8. 

1 

Xp = lim Xp(fcr) = lim — V ln7pi, (79) 

K — >00 K — CXj rvl 

i=l 

where ^ji, j = 1,2, ... ,p, i — 1,2,... are quantities evaluated during the 
successive orthonormahzation procedures (equations ([75]) and ^Bi). Note that 
for p = 1 equation ([79]) is actually equation ([5^ with = 71^. 



6.2 The numerical algorithm for the standard method 

In practice, in order to compute the p largest 1-LCEs with 1 < p < 2N 
we follow the evolution of p initially orthonormal deviation vectors Wj(0) — 
Wj(0) and every t = t time units we replace the evolved vectors ■Wj{kT) 
j = 1, 2, . . . ,p, k = 1,2,... by a new set of orthonormal vectors produced 
by the Gram-Schmidt orthonormahzation method (|76p . During the orthonor- 
mahzation procedure the quantities 'jjk are computed and X1jX27 • ■ • jXp are 
estimated from equation (|79p . This algorithm is described in pseudo-code in 
Table [2] and can be used for the computation of few or even all 1-LCEs. A 
Fortran code of this algorithm can be found in [1441, while [117\ contains 
a similar code developed for the computer algebra platform "Mathematica" 
(Wolfram Research Inc.). 

Let us illustrate the implementation of this algorithm in the particular 
case of the computation of the 2 largest LCEs xi and X2- As shown in Figure 
[8] we start our computation with two orthonormal deviation vectors wi(0) 
and W2(0) which are evolved to wi(t), W2(t) aX t ~ t. Then according 
to the the Gram-Schmidt orthonormahzation method (|75p we define vectors 
Ui(r) and U2(t). In particular, Ui(r) coincides with 'Wi{t) while, U2(r) is the 
component of vector W2(t) in the direction perpendicular to vector Ui(t). 
The norms of these two vectors define the quantities 711 — ||ui(r)||, 721 = 
||u2(t)|| needed for the estimation of xi, X2 from equation ([79|) . Then vectors 
Wi(r) and W2(t) are defined as unitary vectors in the directions of Ui(t) 
and U2(t) respectively. Since the unitary vectors Wi(t), W2(t) are normal by 
construction they constitute the initial set of orthonormal vectors for the next 
iteration of the algorithm. From Figure[8]we easily see that the parallelograms 
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Table 2. The standard method. The algorithm for the computation of the p 
largest LCEs xi , X2 , ■ • • , Xp ''■s limits for t ^ 00 of quantities X\ (t) , X2 (t), . . . , Xp (t) 
l|7ip . according to equation (|79p . The program computes the evolution of 
Xi{t),X2{t), . . . ,Xp{t) with respect to time f up to a given upper value of time 
t = Tm or until any of the quantities Xi{t),X2{t), . . . ,Xp{t) attain a very small 
value, smaller than a low threshold value Xm- 

Input: 1. Hamilton equations of motion ((2]) and variational equations ([8}, or 
equations of the map ^ and of the tangent map (jlip . 

2. Number of desired LCEs p. 

3. Initial condition for the orbit x(0). 

4. Initial orthonormal deviation vectors Wi(0), W2(0), . . ., Wp(0). 

5. Renormalization time r. 

6. Maximal time: Tm and minimum allowed value of Xi{t), 

X2{t), Xpit): Xm. 

Step 1 Set the stopping flag, SF ^ 0, and the counter, <— 1. 

Step 2 While {SF = 0) Do 

Evolve the orbit and the deviation vectors from time t = {k — l)r 
to t = kr, i. e. Compute x{kT) and wi(fcT), W2(fcr), . . ., Wp(fcT). 
Step 3 Perform the Gram-Schmidt orthonormalization procedure 
according to equation (|76l) : 
Do for j = 1 to p 

Compute current vectors Uj(fcr) and values of 7^^. 

Compute and Store current values of Xjikr) = X^jLi ln7ji- 

Set -WjikT) ^ Uj{kT) /"fjk- 
End Do 

Step 4 Set the counter k ^ k + 1. 

Step 5 If [(fcr > Tm) or (Any of Xj((fc - l)r) < X^, j = 1,2, . . . ,p)] Then 
Set SF^ 1. 
End If 
End While 

Step 6 Report the time evolution of Xi (t) , X2 (t), . . . , Xp (t) . 



defined by vectors Wi(r), W2(t) and by vectors Ui(t) and U2(t) have the same 
area. This equality corresponds to the particular case p = 2, i = 1 of equation 
(|77|) . Evidently, since vectors Ui(t), U2(t) are perpendicular to each other, 
we have V0I2 (ui(t), U2(t)) = 711721 in accordance to equation ([75]) . 

6.3 Connection between the standard method and the QR 
decomposition 

Let us rewrite equations (j75p of the Gram-Schmidt orthonormalization proce- 
dure, by solving them with respect to Wj(iT), j = 1, 2, . . . ,p, with 1 < p < 2N 



wi(ir) 
W2(ir) 
W3(ir) 



-7nwi(ir) (80) 
= (wi(iT), W2(ir)}wi(iT) -h72iW2(ir) 

(wi(iT), W3(ir)}wi(iT) -|- (w2(iT), W3(ir)}w2(iT) -f 73iW3(ir) 
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Fig. 8. Numerical sclieme for the computation of the 2 largest LCEs xi, X2 accord- 
ing to the standard method. The orthonormal deviation vectors wi(0), W2(0) are 
evolved according to the variational equations (|SJ (continuous time) or the equations 
of the tangent map (discrete time) for t = t time units. The evolved vectors 
wi (r ) , W2 (r ) , are replaced by a set of orthonormal vectors wi (r ) , W2 (r ) , which span 
the same 2-dimensional vector space, according to the Gram-Schmidt orthonormal- 
ization method (|76p . Then these vectors are again evolved and the same procedure 
is iteratively applied. For each successive time interval [(i — l)r, ir], i = 1, 2, . . ., the 
quantities 71; — ||ui(ir)||, ^21 ~ ||u2(ir)|| are computed and Xi, X2 are estimated 
from equation (|79|l . 



and get the general form 

^k{iT) = ^{wj{iT),Wk{iT))<Vj{iT) + jki<^k{iT), fc = 1,2, . . . ,p. 

This set of equations can be rewritten in matrix form as follows: 

[wi(zr) W2(iT) • • • Wp(zT)] = [wi(zr) W2(«r) • • • Wp(ir)] • 
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7ii (wi(ir), W2(ir)) (wi(ir), W3(ir)} 
72 








(w2(ir), W3(ir)} • ■ 
73i 



(wi(iT),Wp(ir)} 
(w2(iT),Wp(ir)) 
(w3(iT),Wp(ir)) 



7p« 



So the 2iV xp matrix W(ir) — [wi(ir) W2(iT) • • • Wp(iT) ] , having as columns 
the hnearly independent deviation vectors Wj(ir), j = 1, 2, . . . ,p is written 
as a product of the 2N x p matrix Q = [wi(ir) W2(iT) • • ■ Wp(ir)], having 
as columns the coordinates of the orthonormal vectors Wj(ir), j — 1,2, ... ,p 
and satisfying Q^'^Q = Ip, and of an upper triangular pxp matrix R(ir) with 
positive diagonal elements 

= 7r" j '2, i = l,2,.... 

From equations (j80[) we easily see that {w j (it) , w j (ir)) — 'jji and so matrix 
R(zt) can be also expressed as 



R(ir) 



(wi(iT), vifi(iT)) (wi(iT), vif2(iT)) •■■ (wi(ir), wrp(iT)) 
(w2(iT), W2(ir)) •■• (w2(ir), Wp(iT)) 











(Wp(zT),V^fp(iT)) 



The above procedure is the so-called QR decomposition of a matrix. In 
practice, we proved by actually constructing the Q and R matrices via the 
Gram-Schmidt orthonormalization method, the following theorem: 

Theorem 4. Let A be an n x m (n > m) matrix with linearly independent 
columns. Then A can be uniquely factorized as 

A = Q R, 

where Q is an n x m matrix with orthogonal columns, satisfying Q = /,„ 
and R is an m X m invertible upper triangular matrix with positive diagonal 
entries. 

Although we presented the QR decomposition through the Gram-Schmidt 
orthonormalization procedure this decomposition can also be achieved by oth- 
ers, computationally more efficient techniques like for example the House- 
holder transformation [62] [1071 §2.10]. 

Observing that the quantities jji, j — 1,2 ... ,p, i — 1,2..., needed for 
the evaluation of the LCEs through equation (jTO]) are the diagonal elements of 
R(zt) we can implement a variant of the standard method for the computation 
on the LCEs, which is based on the QR decomposition procedure jJH ^21 [SSI 
|40] . Similarly to the procedure followed in Section lOl in order to compute the 
p {1 < p < 2N) largest LCEs we follow the evolution of p initially orthonormal 
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deviation vectors Wj(0) = Wj(0), j ~ 1, 2 . . . ,p, which can be considered as 
cohimns of a 2N x p matrix Q(0). Every t = t time units the matrix W(ir), 
i = 1, 2, . . ., having as columns the deviation vectors 



)'^^wj((z-l)r) =w,(zr), j = l,2, 



"x((i-l)T)'» 

i. e. the columns of Q((i — 1)t) evolved in time interval [(i — l)r, ir] by the 
action of dx((i-i)T)^^i undergoes the QR decomposition procedure 

W(ir) = Q(iT) • TL{iT) (81) 

and the new Q(ir) is again evolved for the next time interval [ir, (i + l)r], 
and so on and so forth. Then the LCEs are estimated from the values of the 
diagonal elements of matrix R(ir) as 



1 

1=1 



(82) 



The corresponding algorithm is presented in pseudo-code in Tabled From the 
above-presented analysis it becomes evident that the standard method devel- 
oped by Shimada and Nagashima |119j and Benettin et al. [14 for the compu- 
tation of the LCEs, is practically a QR decomposition procedure performed by 
the Gram-Schmidt orthonormalization method, although the authors of these 
papers formally do not refer to the QR decomposition. We note that both the 
standard method and the QR decomposition technique presented here can be 
used for the computation of part (p < 2iV) or of the whole {p — 2N) spectrum 
of LCEs. 

As a final remark on the QR decomposition technique let us show the va- 
lidity of equation ((77|) by considering the QR decomposition of matrix W(ir) 
^Tj). According to equations (|105p and (|106p we have 



det W'(iT) -Wiir] 



det (li^iiT) ■ Q^iir) ■ Qiir) ■ R(iT) 

|detR(iT)| 
p 



'detR^(iT) detR(ir) 
^l|u,(^r)|| = 

where the identities Q'^Q = Ip and detR(ir) = 0^=1 7ji have been used. 



p 



6.4 Other methods for computing LCEs 



Over the years several methods have been proposed and applied for computing 
the numerical values of the LCEs. The standard method we discussed so 
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Table 3. Discrete QR decomposition. The algorithm for the computation of 
the p largest LCEs Xii X2, ■ • ■ , Xp according to the QR decomposition method. The 
program computes the evolution of X\[t) , X2{t) , . . . , Xp{t) with respect to time t up 
to a given upper value of time t = Tm or until any of the these quantities becomes 
smaller than a low threshold value Xm- 

Input: 1. Hamilton equations of motion ([2]) and variational equations ((8)1, or 
equations of the map (|4]) and of the tangent map 

2. Number of desired LCEs p. 

3. Initial condition for the orbit x(0). 

4. Initial matrix Q(0) having as columns orthonormal deviation vectors 

Wl(0), W2(0), .., Wp(0). 

5. Time interval r between successive QR decompositions. 

6. Maximal time: Tm and minimum allowed value of Xi{t), 

Xijt), Xpjt): Xm. 

Step 1 Set the stopping flag, SF ^ 0, and the counter, <— 1. 

Step 2 While (SF = 0) Do 

Evolve the orbit and the matrix Q((fc — l)r) from time t = (k — l)r 

to t = fcr, i. e. Compute x.{kT) and W(ir). 
Step 3 Perform the QR decomposition of W(ir) according to (|81|) : 

Compute Q(fcr) and R(fcr). 

Compute and Store current values of Xj(fcr) = X^iLi In R-ij (*''') i 
j = l,2...,p. 
Step 4 Set the counter k ^ fc + 1. 

Step 5 If [(fer > Tm) or (Any of Xj ((fc - l)r) < j = 1, 2, . . . ,p)] Then 
Set SF^ 1. 
End If 
End While 

Step 6 Report the time evolution of Xi (t) , X2 (t) , . . . , Xp (t) . 



far, is the first and probably the simplest method to address this problem. 
As we showed in Section [^31 the standard method, which requires successive 
applications of the Gram-Schmidt orthonormalization procedure, is practically 
equivalent to the QR decomposition technique. 

The reorthonormalization of deviation vectors plays an indispensable role 
for computing the LCEs and the corresponding methods can be distinguished 
in discrete and continuous methods. The discrete methods iteratively approx- 
imate the LCEs in a finite number of (discrete) time steps and therefore 
apply to both continuous and discrete dynamical systems [62l [36l |40] . The 
standard method and its QR decomposition version, are discrete methods. 
A method is called continuous when all relevant quantities are obtained as 
solutions of certain ordinary differential equations, which maintain orthonor- 
mality of deviation vectors continuously. Therefore such methods can only be 
formulated for continuous dynamical systems and not for maps. The use of 
continuous orthonormalization for the numerical computation of LCEs was 
first proposed by Goldhirsch et al. [63] and afterwards developed by several 
authors [67l [62l [Ml 1401 IMl [TTOl [T09l l94l [88] . 
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Discrete and continuous methods are based on appropriate decomposition 
of matrices performed usually by the QR decomposition or by the SVD pro- 
cedure. The discrete QR decomposition, which has already been presented in 
Section 16.31 is the most frequently used method and has proved to be quite 
efficient and reliable. The continuous QR decomposition, as well as methods 
based on the SVD procedure are discussed in some detail at the end of the 
current section. 

Variants of these techniques have been also proposed by several authors. 
Let us briefly refer to some of them. Rangarajan et al. [110] introduced a 
method for the computation of part or of the whole spectrum of LCEs for 
continuous dynamical systems, which does not require rescaling and renor- 
malization of vectors. The key feature of their approach is the use of explicit 
group theoretical representations of orthogonal matrices, which leads to a set 
of coupled ordinary differential equations for the LCEs along with the various 
angles parameterizing the orthogonal matrices involved in the process. Rama- 
subramanian and Sriram (l09j showed that the method is competitive with 
the standard method and the continuous QR decomposition. 

Carbonell et al. [20] proposed a method for the evaluation of the whole 
spectrum of LCEs by approximating the differential equations describing the 
evolution of an orbit of a continuous dynamical system and their associated 
variational equations by two piecewise linear sets of ordinary differential equa- 
tions. Then an SVD or a QR decomposition-based method is applied to these 
two new sets of equations, allowing us to obtain approximations of the LCEs 
of the original system. An advantage of this method is that it does not require 
the simultaneous integration of the two sets of piecewise linear equations. 

Lu et al. [94J proposed a new continuous method for the computation 
of few or of all LCEs, which is related to the QR decomposition technique. 
According to their method one follows the evolution of orthogonal vectors, 
similarly to the QR method, but does not require them to be necessarily 
orthonormal. By relaxing the length requirement Lu et al. [M] established a 
set of recursive differential equations for the evolution of these vectors. Using 
symplectic Runge-Kutta integration schemes for the evolution of these vectors 
they succeeded in preserving automatically the orthogonality between any two 
successive vectors. Normalization of vectors occurs whenever the magnitude 
of any vector exceeds given lower or upper bounds. 

Chen et al. [21] proposed a simple discrete QR algorithm for the computa- 
tion of the whole spectrum of LCEs of a continuous dynamical system. Their 
method is based on a suitable approximation of the solution of variational 
equations by assuming that the Jacobian matrix remains constant over small 
integration time steps. Thus, the scheme requires the numerical solution of the 
2N equations of motion but not the solution of the (2N)'^ variational equa- 
tions since their solution is approximated by an explicit expression involving 
the computed orbit. This approach led to a computationally fast evaluation 
of the LCEs for various multidimensional dynamical systems studied in [24j . 
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It is worth mentioning here a completely different approach, with respect 
to the above-mentioned techniques, which was developed at the early 80's. 
In particular, Fr0yland proposed in [^D] an algorithm for the computation of 
LCEs, which he claimed to be quite efficient in the case of low-dimensional 
systems, and applied it to the Lorenz system [61]. The basic idea behind this 
algorithm is the implementation of appropriate differential equations describ- 
ing the time evolution of volume elements around the orbits of the dynamical 
system, instead of defining these volumes through deviation vectors whose 
evolution is governed by the usual variational equations 

Apart from the actual numerical computation of the values of the LCEs, 
methods for the theoretical estimation of those values have been also devel- 
oped. For example, Li and Chen [9^ provided a theorem for the estimation 
of lower and upper bounds for the values of all LCEs in the case of discrete 
maps. These results were also generalized for the case of continues dynamical 
systems [M]. The validity of these estimates was demonstrated by a compar- 
ison between the estimated bounds and the numerically computed spectrum 
of LCEs of some specific dynamical systems [90l El] . 

Finally, let us refer to a powerful analytical method which allows one to 
verify the existence of positive LCEs for a dynamical system, the so-called 
cone technique. The method was suggested by Wojtkowski [142j and has been 
extensively applied for the study of chaotic billiards |142[ I143[ l43l [97] and 
geodesic flows [41] |42l [19] . A concise description of the techniques can also be 
found in [7] and [25] §3.13]. Considering the space R" a cone C^, with 7 > 0, 
centered around M."~'' is 

= {(u,v) e M'= X M"-*^ : ||u|| < 7||v||} U (0,0). (83) 

Note that {0} x R"-^ c for every 7. In the particular case of n = 3, 
k = 2, C-y corresponds to the usual 3~dimensional cone, while in the case of 
the plane (n = 2) a cone around an axis L is the set of vectors of 
that make angle (f> < arctan7 with the line L. In the case of Hamiltonian 
systems (and symplectic maps) a cone can get the simple form (5q • 5p > 0. 
Finding an invariant family of cones (|83p in T^S, which are mapped strictly 
into themselves by d^'l'* , guarantees that the values of the n — k largest LCEs 
are positive [142) I143j . We emphasize that the cone technique is not used for 
the explicit numerical computation of the LCEs, but for the analytical proof 
of the existence of positive LCEs, providing at the same time some bounds 
for their actual values. 

Continuous QR decomposition methods 

The QR decomposition methods allow the computation of all or of the p 
{1 < p < 2N) largest LCEs. Let us discuss in more detail the developed 
procedure for both cases following mainly [62] [36] [94] . 
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Computing the complete spectrum of LCEs 

The basic idea of the method is to avoid directly solving the differential equa- 
tion pO)) . by requiring Y{t) ~ Q(t)R(t) where Q{t) is orthogonal and R(t) 
is upper triangular with positive diagonal elements, according to Theorem 3) 
With this decomposition, one can write equation (jlOp into the form 

Q^Q + RR^ = Q^AQ, 

where, for convenience, we dropped out the explicit dependence of the matrices 
on time t, i. e. Q{t) = Q. Since Q^'^Q is skew and RR^^ is upper triangular, 
one reads off the differential equations 

Q = QS, (84) 

where S is the skew symmetric matrix 

S =^ Q^Q 

with elements 

r (Q^AQ),, i>j 
Sy = <^ i^j , i,j = l,2,...,2iV, (85) 

[-(Q^AQ),, i<j 

and 

1^ - (Q^AQ)pp, p, - 1, 2, . . . , 27V (86) 

XVpp 

where Rpp are the diagonal elements of R. As we have already seen in equation 
(|82|l the LCEs are related to the elements Rpp, through 

Xp = lim i In Rpp (i). 

t^oo t 

Thus, in order to compute the spectrum of LCEs only equations ([M)) and ([5^ 
have to be solved simultaneously with the equations of motion In practice, 
the knowledge of matrix R is not necessary for the actual computation of the 
LCEs. Noticing that 

-^(lnRpp) = ^ = (QTAQ)pp = qp.Aqp, (87) 

where qp is the pth column vector of Q, we can compute the LCEs using 

1 /•* 

Xp = lim - / qp ■ Aqpdt. 

t-*ca t Jq 

In practice, the LCEs can be estimated through a recursive formula. Let 
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Then we have 



1 

KT Jq 



dt. 



1 



Xp {{k + l)r) = J^^-Jy qp • Mpdt 

I fkr ^ /■{k+l)T 



(fc+1) 



Replacing the first integral with kTXp{kT) we get 



Xp {(k + l)r) - ^Xpikr) + j^^-jy^ ^ qp • Aq^dt, (88) 



and 



= lim Xp{kT). (89) 

/c — ^oc 

The basic difference between the discrete QR decomposition method pre- 
sented in Section [^751 and the continuous QR method presented here, is that 
in the first method the orthonormalization is performed numerically at dis- 
crete time steps, while the latter method seeks to maintain the orthogonality 
via solving differential equations that encode the orthogonality continously. 

Computation of the p > 1 largest LCEs 

If we want to compute the p largest LCEs, with 1 < p < 2A'', we change 
equation (fTO|) to 

Y{t) ^ A{t)Y{t) , with Y(0)'^Y(0) = Ip. (90) 

where Y(t) is in practice, the 2N x p matrix having as columns the p deviation 
vectors wi{t),W2{t), . . . ,Wp{t). Applying Theorem|4]we get Y{t) = Q(t)R(i) 
where Q{t) is orthogonal so that the identity Q^'^Q = I holds but not the 
QQ"^ = I. Then from equation ([90|) we get 

R = (q'^AQ - s) R 

where S = Q'^Q is a p x p matrix whose elements are given by equation 
for i, j — 1,2, . . . ,p. Since R is invertible, from the relations 



and 

we obtain 



RR^ = Q^AQ S 



Q = AQ QRR \ (91) 
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Q = (a QQ^A + QSQt) Q, 

or 

Q = H(Q,i)Q, (92) 

with 

H(Q, t)^A- QQ^A + QSQ^. 

Notice that the matrix H(Q,t) in not necessarily skew-symmetric, and 
the term QQ"^ is responsible for lack of skew-symmetry in H. Of course for 
p = 2N equation ([9^ reduces to equation Q = QS (IMl) . The evolution of the 
diagonal elements of R are again governed by equation ([55]) . but for p < 2N, 
and so the p largest LCEs can be computed again from equations (|87p - (|89|) . 

The main difference with respect to the case of the computation of the 
whole spectrum is the numerical difficulties arising in solving equation (|92p . 
since H is not skew-symmetric as was matrix S in equation (|84p . Due to this 
difference usual numerical integration techniques fail to preserve the orthog- 
onality of matrix Q. 

A central observation of [36 is that the matrix H has a weak skew- 
symmetry property. The matrix H is called weak skew-symmetric if 

Q'^ (h(Q, t) + H'^(Q, i)) Q = 0, whenever Q^Q = Ip. 

A matrix H is said to be strongly skew-symmetric if it is skew-symmetric, 
i. e. H""" — — H. Christiansen and Rugh [26j proposed a method according to 
which, the numerically unstable equations (PT|) for the continuous orthonor- 
malization could be stabilized by the addition of an appropriate dissipation 
term. This idea was also used in [18] , where it was shown that it is possible to 
reformulate equation (I92|l so that H becomes strongly skew-symmetric and 
thus, achieve a numerically stable algorithm for the computation of few LCEs. 

Discrete and continuous methods based on the SVD procedure 

An alternative way of evaluating the LCEs is obtained by applying the SVD 
procedure on the fundamental 2N x 2N matrix Y(t), which defines the evo- 
lution of deviation vectors through equations ([9]) and p2|) for continuous and 
discrete systems respectively. According to the SVD algorithm a 2Nxp matrix 
{p < 2N) B can be written as the product of a 2N xp column-orthogonal ma- 
trix U, apxp diagonal matrix F with positive or zero elements ai, i — 1, . . . ,p 
(the so-called singular values), and the transpose of a, pxp orthogonal matrix 
V: 

B = U F V^. 

We note that matrices U and V are orthogonal so that: 



U'^ ■ U = V'^ ■ V = L 



(93) 
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For a more detailed description of the SVD method, as well as an algorithm for 
its implementation the reader is referred to [107[ Section 2.6] and references 
therein. The SVD is unique up to permutations of corresponding columns, 
rows and diagonal elements of matrices U, V and F. Advanced numerical 
techniques for the computation of the singular values of a product of many 
matrices can be found for example in [1301 llOlj . 
So, for the purposes of our study let 

Y U ■ F ■ V"^, (94) 

where we dropped out as before, the explicit dependence of the matrices on 
time t. In those cases where all singular values are different, a unique decom- 
position can be achieved by the additional request of a strictly monotonically 
decreasing singular value spectrum, i. e. <7i{t) > a2{t) > ■ ■ ■ > (J2N{i)- Multi- 
plying equation (|94|) with the transpose 

= V ■ F'^ • U'^, 

from the left we get 

• Y ^ V • F^ • U"^ • U ■ F ■ ^ V • diag(cr,f (t)) • V^, (95) 

where equation ([55]) has been used. From equation (|95p we see that the eigen- 
values of the diagonal matrix diag(crf (i)), i. e. the squares of the singular 
values of Y(<), are equal to the eigenvalues of the symmetric matrix Y'^Y. 
Then from point 4 of the MET we conclude that the LCEs are related to the 
singular values of Y{t) through [62l I130j 

Xp - lim ilna,(t), p = 1, 2, . . . , 2iV, 

which implies that the LCEs can be evaluated as the limits for i ^ cx3 of the 
time rate of the logarithms of the singular values. 

Theoretical aspects of the SVD technique, as well as a detailed study of 
its ability to approximate the spectrum of LCEs can be found in f lOll [571 [55] . 
Continuous [67l [62l [39] and discrete [130] versions of the SVD algorithm have 
been applied for the computation of few or of all LCEs, although this approach 
is not widely used. A basic problem of these methods is that they fail to 
compute the spectrum of LCEs if it is degenerate, i. e. when two or more 
LCEs are equal or very close to each other, due to the appearance of ill- 
conditioned matrices. 



7 Chaos detection techniques 

A simple, qualitative way of studying the dynamics of a Hamiltonian system 
is by plotting the successive intersections of its orbits with a Poincare surface 
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of section (PSS) (e. g. [72] [HH p. 17-20]). Similarly, in the case of symplectic 
maps one simply plots the phase space of the system. This qualitative method 
has been extensively applied to 2d maps and to 2D Hamiltonians, since in 
these systems the PSS is a 2-dimensional plane. In such systems one can 
visually distinguish between regular and chaotic orbits since the points of a 
regular orbit lie on a torus and form a smooth closed curve, while the points 
of a chaotic orbit appear randomly scattered. In 3D Hamiltonian systems (or 
4d symplectic maps) however, the PSS (or the phase space) is 4-dimensional 
and the behavior of the orbits cannot be easily visualized. Things become 
even more difficult and deceiving for multidimensional systems. One way to 
overcome this problem is to project the PSS (or the phase space) to spaces 
with lower dimensions (see e.g. |139| 11401 I105j ) although these projections 
are often very complicated and difficult to interpret. Thus, we need fast and 
accurate numerical tools to give us information about the regular or chaotic 
character of orbits, mainly when the dynamical system has many degrees of 
freedom. 

The most commonly employed method for distinguishing between regular 
and chaotic behavior is the evaluation of the mLCE xi, because if xi > the 
orbit is chaotic. The main problem of using the value of xi as an indicator 
of chaoticity is that, in practice, the numerical computation may take a huge 
amount of time, in particular for orbits which stick to regular ones for a 
long time before showing their chaotic behavior. Since xi is defined as the 
limit for i — > cxi of the quantity Xi{t) (j54|) . the time needed for Xi{t) to 
converge to its limiting value is not known a priori and may become extremely 
long. Nevertheless, we should keep in mind that the mLCE gives us more 
information than just characterizing an orbit as regular or chaotic, since it 
also quantifies the notion of chaoticity by providing a characteristic time scale 
for the studied dynamical system, namely the Lyapunov time (|51|) . 

In order to address the problem of the fast and reliable determination of 
the regular or chaotic nature of orbits, several methods have been developed 
over the years with varying degrees of success. These methods can be divided 
in two major categories: Some are based on the study of the evolution of 
deviation vectors from a given orbit, like the computation of xi, while others 
rely on the analysis of the particular orbit itself. 

Among other chaoticity detectors, belonging to the same category with 
the evaluation of the mLCE, are the fast Lyapunov indicator (FLI) [5S1 |M1 
[56l [89l l49l [69] and its variants [H [5], the smaller alignment index (SALI) 
|1221 I124|, I125j and its generalization, the so-called generalized alignment 
index (GALI) [1261 I127j . the mean exponential growth of nearby orbits 
(MEGNO) [211 [29], the relative Lyapunov indicator (RLI) [113 [HI], the aver- 
age power law exponent (APLE) [SS], as well as methods based on the study 
of spectra of quantities related to the deviation vectors like the stretching 
numbers [ST] [93] 11351 1138] , the helicity angles (the angles of deviation vectors 
with a fixed direction) [32], the twist angles (the differences of two succes- 
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sive helicity angles) [35] , or the study of the differences between such spectra 



In the category of methods based on the analysis of a time series con- 
structed by the coordinates of the orbit under study, one may list the frequency 
map analysis of Laskar [83l [M [Ml IS , the '0-1' test [Ml El], the method of 
the low frequency spectral analysis |137l[8T] . the 'patterns method' [1201ll21j . 
the recurrence plots technique [1471 1148j and the information entropy index 
|100j . One could also refer to several ideas presented by various authors that 
could be used in order to distinguish between chaoticity and regularity, like 
the differences appearing for regular and chaotic orbits in the time evolutions 
of their correlation dimension [50j . in the time averages of kinetic energies 
related to the virial theorem [74] and in the statistical properties of the series 
of time intervals between successive intersections of orbits with a PSS [80] , 

A systematic and detailed comparative study of the efficiency and reliabil- 
ity of the various chaos detection techniques has not been done yet, although 
comparisons between some of the existing methods have been performed spo- 



radically in studies of particular dynamical systems [122\ I125[ 11321 11331 l82l 

[MIE]. 



Let us now focus our attention on the behavior of the FLI and of the GALI 
and on their connection to the LCEs. The FLI was introduced as an indicator 
of chaos in [58l[59j and after some minor modifications in its definition, it was 
used for the distinction between resonant and not resonant regular motion 
[561149]. The FLI is defined as 



where vf(t) is a deviation vector from the studied orbit at point x(i), which 
initially had unit norm, i. e. ||w(0)|| = 1. In practice, FLI(i) registers the 
maximum length that an initially unitary deviation vector attains from the 
beginning of its evolution up to the current time t. Using the notation ap- 
pearing in equation (|59|) . the FLI can be computed as 



with the initial norm Dq of the deviation vector being Dq = 1. 

According to equation (|62p the norm of w(i) increases linearly in time in 
the case of regular orbits. On the other hand, in the case of chaotic orbits the 
norm of any deviation vector exhibits an exponential increase in time, with 
an exponent which approximates xi for t —>■ oo. Thus, the norm of a deviation 
vector reaches rapidly completely different values for regular and chaotic or- 
bits, which actually differ by many orders of magnitude. This behavior allows 
FLI to discriminate between regular orbits, for which FLI has relatively small 
values, and chaotic orbits, for which FLI gets very large values. 

The main difference of FLI with respect to the evaluation of the mLCE 
by equation (|59p is that FLI registers the current value of the norm of the 



FLI(i) = sup ln|lw(t)||. 
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deviation vector and does not try to compute the limit value, for t —>■ oo, 
of the mean of stretching numbers as xi does. By dropping the time average 
requirement of the stretching numbers, FLI succeeds in determining the nature 
of orbits faster than the computation of the mLCE. 

The generalized alignment index of order p (GALIp) is determined through 
the evolution of 2 < p < 2N initially linearly independent deviation vectors 
Wi(0), i = 1,2, ... ,p and so it is more related to the computation of many 
LCEs than to the computation of the mLCE. The evolved deviation vectors 
MVilt) are normalized from time to time in order to avoid overflow problems, 
but their directions are left intact. Then, according to 126 GALIp is defined to 
be the volume of the p-parallelogram having as edges the p unitary deviation 
vectors Wi^t), i — 1,2, ... ,p 

GALIp(i) = ||wi(t) Aw2(t) A--- Awp(i)||. (96) 

In [126] the value of GALIp is computed according to equation (|105p . while 
in [21 1127j a more cfhcient numerical technique based on the SVD algorithm 
is applied. From the definition of GALIp it becomes evident that if at least 
two of the deviation vectors become linearly dependent, the wedge product in 
(|M)) becomes zero and the GALIp vanishes. 

In the case of a chaotic orbit all deviation vectors tend to become linearly 
dependent, aligning in the direction which corresponds to the mLCE and 
GALIp tends to zero exponentially following the law |126j : 

GALIp(i) - e"[(xi-X2)+(xi-X3)+-+(xi-Xp)]*^ 

where Xi: X2: ■ ■ • j Xp ^-re the p largest LCEs. On the other hand, in the case 
of regular motion all deviation vectors tend to fall on the iV-dimensional 
tangent space of the torus on which the motion lies. Thus, if we start with 
p < N general deviation vectors they will remain linearly independent on the 
iV-dimensional tangent space of the torus, since there is no particular reason 
for them to become linearly dependent. As a consequence GALIp remains 
practically constant for p < N. On the other hand, GALIp tends to zero for p > 
N, since some deviation vectors will eventually become linearly dependent, 
following a particular power law which depends on the dimensionality N of 
the torus and the number p of deviation vectors. So, the generic behavior of 
GALIp for regular orbits lying on A^-dimensional tori is given by |126j 

^»xT /,\ f constant if 2 < p < 

^^^^^(^)-|^ ilN-<p-<2N- (97) 

The different behavior of GALIp for regular orbits, where it remains dif- 
ferent from zero or tends to zero following a power law, and for chaotic or- 
bits, where it tends exponentially to zero, makes GALIp an ideal indicator of 
chaoticity independent of the dimensions of the system [1261 11271 [T5] . GALIp 
is a generalization of the SALI method [1221 11241 1125j which is related to the 
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evolution of only two deviation vectors. Actually GALI2 oc SALI. However, 
GALIp provides significantly more detailed information on the local dynamics, 
and allows for a faster and clearer distinction between order and chaos. It was 
shown recently [571 11^7] that GALIj, can also be used for the determination 
of the dimensionality of the torus on which regular motion occurs. 

As we discussed in Section EH] the alignment of all deviation vectors to the 
direction corresponding to the niLCE is a basic problem for the computation 
of many LCEs, which is overcome by successive orthonormalizations of the 
set of deviation vectors. The GALIs on the other hand, exploit exactly this 
'problem' in order to determine rapidly and with certainty the regular or 
chaotic nature of orbits. 

It was shown in Section 14.11 that the values of all LCEs (and therefore 
the value of the mLCE) do not depend on the particular used norm. On the 
other hand, the quantitative results of all chaos detection techniques based 
on quantities related to the dynamics of the tangent space on a finite time, 
depend on the used norm, or on the coordinates of the studied system. For 
example, the actual values of the finite time mLCE Xi{t) (jM]) will be different 
for different norms, or for different coordinates, although its limiting value for 
t ^ 00, i. e. the mLCE xi, will be always the same. Other chaos detection 
methods, like the FLI and the GALI, which depend on the current values of 
some norm-related quantities and not on their limiting values for t ^ 00, 
will attain different values for different norms and/or coordinate systems. 
Although the values of these indices will be different, one could expect that 
their qualitative behavior would be independent of the chosen norm and the 
used coordinates, since these indices depend on the geometrical properties of 
the deviation vectors. For example, the GALI quantifies the linear dependence 
or independence of deviation vectors, a property which obviously does not 
depend on the particular used norm or coordinates. Indeed, some arguments 
explaining the independence of the behavior of the GALI method on the 
chosen coordinates can be found in [126j . Nevertheless, a systematic study 
focused on the influence of the used norm on the qualitative behavior of the 
various chaos indicators has not been performed yet, although it would be of 
great interest. 

8 LCEs of dissipative systems and time series 

The presentation of the LCEs in this report was mainly done in connection 
to conservative dynamical systems, i. e. autonomous Hamiltonian flows and 
symplectic maps. The restriction to conservative systems is not necessary since 
the theory of LCEs, as well as the techniques for their evaluation are valid 
for general dynamical systems like for example dissipative ones. In addition, 
within what is called time series analysis (see e.g. [TSj) it is of great interest to 
measure LCEs in order to understand the underlying dynamics that produces 
any time series of experimental data. For the completeness of our presentation 
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we devote the last section of our report to a concise survey of results concerning 
the LCEs of dissipative systems and time series. 

8.1 Dissipative systems 

In contrast to Hamiltonian systems and symplectic maps for which the con- 
servation of the phase space volume is a fundamental constraint of the motion, 
a dissipative system is characterized by a decrease of the phase space volume 
with increasing time. This leads to the contraction of motion on a surface of 
lower dimensionality than the original phase space, which is called attractor. 
Thus any dissipative dynamical system will have at least one negative LCE, 
the sum of all its LCEs (which actually measures the contraction rate of the 
phase space volume through equation (|43p) is negative and after some initial 
transient time the motion occurs on an attractor. 

Any continuous n-dimensional dissipative dynamical system without a 
stationary point (which is often called a fixed point) has at least one LCE 
equal to zero [70] as we have already discussed in Section 14.51 For regular 
motion the attractor of dissipative flows represents a fixed point having all 
its LCEs negative, or a quasiperiodic orbit lying on a p-dimensional torus 
(p < n) having p zero LCEs while the rest n — p exponents are negative. For 
dissipative flows in three or more dimensions there can also exist attractors 
having a very complicated geometrical structure which are called 'strange'. 

Strange attractors have one or more positive LCEs implying that the mo- 
tion on them is chaotic. The exponential expansion indicated by a positive 
LCE is incompatible with motion on a bounded attractor unless some sort of 
folding process merges separated orbits. Each positive exponent corresponds 
to a direction in which the system experiences the repeated stretching and 
folding that decorrelates nearby orbits on the attractor. A simple geometri- 
cal construction of a hypothetical strange attractor where orbits are bounded 
despite the fact that nearby orbits diverge exponentially can be found in [92l 
Sect. 1.5]. 

The numerical methods for the evaluation of the mLCE, of the p (1 < p < 
n) largest LCEs and of the whole spectrum of them, presented in Sections 
[S] and [H can be applied also to dissipative systems. Actually, many of these 
techniques were initially used in studies of dissipative models [511 Uli^l [5T1 . 
For a detailed description of the dynamical features of dissipative systems, as 
well as of the behavior of LCEs for such systems the reader is referred, for 
example, to |103[ [44] [92l Sect. 1.5, Chapt. 7 and 8] and references therein. 

8.2 Computing LCEs from a time series 

A basic task in real physical experiments is the understanding of the dynamical 
properties of the studied system by the analysis of some observed time series of 
data. The knowledge of the LCEs of the system is one important step towards 
the fulfillment of this goal. Usually, we have no knowledge of the nonlinear 
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equations that govern the time evolution of the system which produces the 
experimental data. This lack of information makes the computation of the 
spectrum of LCEs of the system a hard and challenging task. 

The methods developed for the determination of the LCEs from a scalar 
time series have as starting point the technique of phase space reconstruction 
with delay coordinates |1041 1134| I112j [78l Chapt. 3 and 9]. This technique is 
used for recreating a d-dimensional phase space to capture the behavior of 
the dynamical system which produces the observed scalar time series. 

Assume that we have Nd measurements of a dynamical quantity x taken 
at times i„ = + '^•''"1 i- e. x{n) = a:(<o + nr), n = 0, 1, 2, . . . , Nd — 1. Then 
we produce Nd = No — (rf — 1)T d-dimensional vectors x(i„) from the a;'s as 



where T is the (integer) delay time. With this procedure we construct Nd 
points in a d-dimensional phase space, which can be treated as successive 
points of a hypothetical orbit. We assume that the evolution of x(t„) to 
x(i„+i) is given by some map and we seek to evaluate the LCEs of this orbit. 

The first algorithm to compute LCEs for a time series was introduced 
by Wolf et al. [144j . According to their method (which is also referred as 
the direct method), in order to compute the mLCE we first locate the nearest 
neighbor (in the Euclidean sense) x(t/j), to the initial point x(to) and define the 
corresponding deviation vector w(to) = x(to) — x(tfc) and its length L{to) = 
||w(to)||- The points x(to) and x(tfc) are considered as initial conditions of 
two nearby orbits and are followed in time. Then the mLCE is evaluated by 
the method discussed in Section [^21 which approximates deviation vectors by 
differences of nearby orbits. So, at some later time (which is fixed a priori 
or determined by some predefined threshold violation of the vector's length) 
the evolved deviation vector vif' (im J = yi{tmi)—^{tk+mi) is normahzed and its 
length L'{tmi) — is registered. The 'normalization' of the evolved 

deviation vector is done by looking for a new data point, say x{ti), whose 
distance L{tmi) = ll^(^mi) ^ x;(^;)|| from the studied orbit is small and the 
corresponding deviation vector w(tm ^ ) — x{t„i^ ) — x(l) has the same direction 
with w'{tmi)- Of course with finite amount of data, one cannot hope to find 
a replacement point x(/) which falls exactly on the direction of w'(t^J but 
chooses a point that comes as close as possible. Assuming that such point is 
found the procedure is repeated and an estimation Xi{tm„) of the mLCE xi 
is obtained by an equation analogous to equation ([56]): 



with mo = 0. A Fortran code of this algorithm with fixed time steps between 
replacements of deviation vectors is given in [144] . 

Generalizing this technique by evolving simultaneously p > 1 deviation 
vectors, i. e. following the evolution of the orbit under study, as well as of p 



x(t„) = [ x{n) x{n + T) . . . 



x{n + {d - l)T)f , 
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nearby orbits, we can, in principle, evaluate the p-mLCE Xi of the system, 
which is equal to the sum of the p largest 1-LCEs (see equation (|57)) V Then 
the values of Xi * = li2, . . . ,p can be computed from equation (|7^ . This 
procedure corresponds to a variant of the standard method for computing 
the LCEs, presented in |119j and discussed in Section [Ql in that deviation 
vectors are defined as differences of neighboring orbits. The implementation 
of this approach requires the repeated replacement of the deviation vectors, 
i. e. the replacement of the p points close to the evolved orbit under consid- 
eration, when the lengths of the vectors exceed some threshold value. This 
replacement should be done in a way that the volume of the corresponding 
p-parallelogram is small, and in particular smaller than the replaced volume, 
and the new p vectors point more or less to the same direction like the old 
ones. This procedure is explained in detail in [144] for the particular case of 

the computation of x]^ — xi + X2, where a triplet of points is involved. 

It is clear that in order to achieve a good replacement of the evolved p 
vectors, which will lead to a reliable estimation of the LCEs, the numerical 
data have to satisfy many conditions. Usually this is not feasible due to the 
limited number of data points. So the direct method of 144J does not yield 
very precise results for the LCEs. Another limitation of the method, which was 
pointed out in Wolf ct al. [144] , is that it should not be used for finding nega- 
tive LCEs which correspond to shrinking directions, due to a cut off in small 
distances implied mainly by the level of noise of the experimental data. An 
additional disadvantage of the direct method is that many parameters which 
influence the estimated values of the LCEs like the embedding dimension rf, 
the delay time T, the tolerances in direction angles during vector replacements 
and the evolution times between replacements, have to be tuned properly in 
order to obtain reliable results. 

A different approach for the computation of the whole spectrum of LCEs is 
based on the numerical determination of matrix Y„, ii — 1,2,..., of equation 
p^ . which defines the evolution of deviation vectors in the reconstructed 
phase space. This method was introduced in [118] and was studied in more 
detail in [44l|45] (see also [78l Chapt. 11]). According to this approach, often 
called the tangent space method^ matrix Y„ is evaluated for each point of 
the studied orbit through local linear fits of the data. In particular, for every 
point x(t„) of the orbit we find all its neighboring points, i. e. points whose 
distance from x(t„) is less than a predefined small value e. Each of these 
point define a deviation vector. Then we find the next iteration of all these 
points and see how these vectors evolve. Keeping only the evolved vectors 
having length less than e we evaluate matrix Y„ through a least-square-error 
algorithm. By repeating this procedure for the whole length of the studied 
orbit we are able to evaluate at each point of the orbit matrix Y„ which defines 
the evolution of deviation vectors over one time step. Then by applying the 
QR decomposition version of the standard method, which was presented in 
Section [6?3l we estimate the values of the LCEs. The corresponding algorithm 
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is included in the TISEAN software package of nonlinear time series analysis 
methods developed by Hegger et al. [7T]. It is also worth mentioning that 
Brown et al. [17j improved the tangent space method by using higher order 
polynomials for the local fit. 

If, on the other hand, we are interested only in the evaluation of the mLCE 
of a time series we can apply the algorithm proposed by Rosenstein et al. Ill] 
and Kantz [77] . The method is based on the statistical study of the evolution of 
distances of neighboring orbits. This approach is in the same spirit of Wolf et 
al. [144] although being simpler since it compares distances and not directions. 
A basic difference with the direct method is that for each point of the reference 
orbit not one, but several neighboring orbits are evaluated leading to improved 
estimates of the mLCE with smaller statistical fluctuations even in the case of 
small data sets. This algorithm is also included in the TISEAN package |71) . 
while its Fortran and C codes can be found in [751 Appendix B]. 
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Appendix 

A Exterior algebra and wedge product: Some basic 
notions 

We present here some basic results of the exterior algebra theory along with 
an introduction to the theory of wedge products following [T] and textbooks 
such as [128[ [68l I129j . We also provide some simple illustrative examples of 
these results. 

Let us consider an M-dimensional vector space V over the field of real 
numbers M. The exterior algebra of V is denoted by A{V) and its multiplica- 
tion, known as the wedge product or the exterior product, is written as A. The 
wedge product is associative: 

(u A v) A w u A (v A w?) 

for u, V, w g y and bilinear 
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(ciU + C2V) A w = Ci(u A w) + C2(v A w), 
w A (ciU + C2V) = ci(w A u) + C2(w A v), 

for u, v,w g y and ci,C2 G M. The wedge product is also alternating on V 

u A u = 

for all vectors u G V. Thus we have that 

u A V = —V A u 

for all vectors u, v e F and 

ui A U2 A • • • A Ufe = (98) 

whenever Ui,U2,...,Ufc e V are linearly dependent. Elements of the form 
Ui A U2 A • ■ • A Ufe with Ui, U2, . . . , Ufc G y are called k-vectors. The subspace 
of A{V) generated by all fc-vectors is called the k-th exterior power of V and 
denoted by A''{V). 

Let {ei, §2, . . . , ei\i} be an orthonormal basis of V, i. e. e^, i = 1, 2, . . . , M 
are linearly independent vectors of unit magnitude and 

where ' • ' denotes the inner product in V and 

„ J 1 for i = j 
" \ for z / J ■ 

It can be easily seen that the set 

{e^^ A A ■ ■ ■ A e.i^ \ 1 < h < 12 < ■ ■ ■ < ik < M} (99) 

is a basis of A^{V) since any wedge product of the form Ui A U2 A • • • A 
can be written as a linear combination of the fc-vectors of equation ([M)) . This 
is true because every vector u^, i = l,2,...,k can be written as a hnear 
combination of the basis vectors e^, i = 1,2, ...,M and using the bilinearity 
of the wedge product this can be expanded to a linear combination of wedge 
products of those basis vectors. Any wedge product in which the same basis 
vector appears more than once is zero, while any wedge product in which the 
basis vectors do not appear in the proper order can be reordered, changing the 
sign whenever two basis vectors change places. The dimension dk oi A''{V) is 
equal to the binomial coefficient 



dk^dimA'^iV) 



M \ Ml 



k J k\{M~-k)\ 



64 Ch. Skokos 



Ordering the elements of basis 
lexicographical order 



of A'^{V) according to the standard 



I < ii < 12 < ■ ■ ■ < ik < M, i 



1,2,- 



any /c-vector u G A''{V) can be represented as 



j=i 



,4, (100) 



(101) 



A fc-vector which can be written as the wedge product of k Hnear independent 
vectors of V is called decomposable. Of course, if the k vectors are linearly 
dependent we get the zero /c-vector (|98p . Note that not all A:-vectors are 
decomposable. For example the 2-vector u = ei A 62 + 63 A 64 € il^(R'*) is 
not decomposable as it cannot be written as Ui A U2 with Ui, U2 G M*. 

Let us consider a decomposable fc-vector u = Ui A U2 A • • • A . Then 
the coefficients Ui in (jlOip are the minors of matrix U having as columns the 
coordinates of vectors Uj, i = 1, 2, . . . , fc with respect to the orthonormal basis 
ei, j = 1, 2, . . . , M. In matrix form we have 



Ui U2 



Ufe = ei 62 



Uii Ui2 

U21 U22 

UMl UM2 



Ulk 
U2k 

UMk 



[ei §2 • • • eA/]-U 



where My, i 



1,2, 



(102) 

, M , j — 1,2, ... ,k are real numbers. Then, the wedge 



product Ui A U2 A • ■ • A Ufc is written as 
u = Ui A U2 A 



E 

l<il<j2<-<ifc<A/ 



A Ufc = ^ UiUJi 

1=1 

Uiil Ui^2 ■ ■ ■ Miifc ' 
Uial Ui22 ■ ■ ■ Ui^k 



(103) 



Ae, 



where the sum is performed over all possible combinations of k indices out of 
the M total indices and | | denotes the determinant. So, the coefficient of a 
particular /c-vector ej^ Ae^^ A - ■ -Ae^^ is the determinant of the kxk submatrix 
of the M X k matrix of coefficients appearing in equation (|102ll formed by its 
ii, 12, ik rows. 

The inner product on V induces an inner product on each vector space 
A'^CV) as follows: Considering two decomposable fc-vectors 



u = Ui A U2 A • • • A Ufe and v = vi A V2 A • • • A v^. 
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with Ui, Vj ^ V, i,j = 1,2, k, the inner product of u, v e A'^{V) is defined 
by 

' Ui ■ Vi Ui ■ V2 • • • Ui ■ Vfc 
U2 • Vi U2 • V2 • • • U2 ■ Vfc 



(u,v), 



dcf 



Ufc • Vi Ufc • V2 • • • Ufc • Vfc 



• V 



(104) 



where U, V are matrices having as columns the coefficients of vectors Uj, Vj, 
i = 1,2, . . . ,k with respect to the orthonormal {ei , §2, . . . , e^} (see equation 
(|102p ). Since every element of A''{V) is a sum of decomposable element, this 
definition extends by bilinearity to any fc-vector. Obviously for the basis (|100[) 
of A''{V) we have 

{u}i,u}.j)k = Sij , i.j = 1,2,.. .,dk. 



implying that the basis is orthonormal. Inner product (|104p defines a norm 
II II for fc- vectors by 



,|u|| = \/(u,u)fe = 



U'^ -U 



Thus, the norm of a decomposable fc-vector ()103|) is given by 

/ W.. X 1/2 



Ui A U2 A • • • A Ufc 



E 

l<ii<i2<-<ik<M 



■ u 



"ill Mii2 • ■ ■ U,^k 



E^? 

2. 1/2 



(105) 



and it measures the volume vol(Pfe) of the /c-parallelogram Pk having as edges 
the k vectors Ui, U2, ■ • • , Ufc, since this volume is defined as (see e. g. [75l 
p. 472]) 



VOl(Pfe) = 



-U 



(106) 



The use of a different orthonormal basis does not change the numerical 
value of vol(Pfc). This can be easily seen as follows: Let f j, i = 1, 2, • • • , M be 
a different orthonormal basis of V related to basis through 

[§1 §2 • • • Gm] = [fi f2 • • ■ hi] ■ A 

where A is an orthogonal M x M matrix, i. e. A^^ = A"*". From equation 
pl}^ we get 

[Ui U2 • • • Ufc] = [fi £2 ■ • • fAf] • A • U, 
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whence the volume vol'(Pfc) with respect to the new basis f^, i = 1,2, 
is given by 

VOl'(Pfc) = 



(A ■ U)^ • A • U 



U"^ • A^i • A • U 



U'^ ■ U 



VOl(Pfc), 



where the orthogonahty of A was used. This resuh is not surprising since an 
orthogonal matrix corresponds to a rotation that leaves unchanged the norms 
of vectors and the angles between them. 
Finally we note that the equality 



E 

l<ii<i2<-<ik<M 



Uiik 



appearing in equation (jl05[) is the so-called Lagrange identity (e. g. |68[ 
p. 108], [M p. 103]). 



A.l An illustrative example 

In order to illustrate the content of the previous section we consider here 
a specific example. Let V be the vector space of M = 4-dimensional real 
vectors, i. e. y = and 

§1 ^ (1, 0, 0, 0) , 62 = (0, 1, 0, 0) , 63 = (0, 0, 1, 0) , 64 = (0, 0, 0, 1) , (107) 

the usual orthonormal basis of R"*. Then the lexicographically ordered or- 
thonormal basis (jlOOp of the 3,2 = 6-dimcnsional vector space yl^(R'') is 



u;i = 61 A 62 , U)2 
u;4 = 62 A 63 , u;5 



§1 A 63 , u;3 = §1 A 64 
§2 A 64 , = 63 A 64 



(108) 



The yl'^(R^) vector space has dimension — A and the set 

yi = 61 A 62 A 63 , y2 = 61 A §2 A 64 , 
y3 = §1 A 63 A 64 , y4 = 62 A 63 A 64 , 

as an orthonormal basis, while the d^ — 1-dimensional vector space yl''(R^) 
is spanned by vector 

xi = 61 A 62 A 63 A 64. 

Let us now consider 4 linearly independent vectors Ui, U2, U3, U4 of R^ 
and the matrix 



U 



Ui U2 U3 U4 



Mil U12 Ui3 Uii 
U2I U22 U23 U24 
U3I W32 M33 W34 
M4I W42 U43 "44 



»,j = 1,2,3,4, 
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having as columns the coordinates of these vectors with respect to the basis 
(fTUTD of Ml 

Considering basis (|108p of yl-^(M^) the 2~vector Ui A U2 is given by 



Ui A U2 = 



Ull 


Ul2 


Wi + 


Ull 


Ul2 




Ull 


"12 










^2 + 






U21 


"22 




U3I 


U32 




U41 


"42 



"21 "22 




"21 "22 




"31 "32 




0)4 + 




^5 + 




"31 "32 




"41 "42 




"41 "42 



^^3 + 



according to equation p03[) . For the norm of this vector we get from equations 
(fT04l) and (fTOS]) : 



I Ui A U2 h 



||Ui||" Ui • U2 
U2 • Ui IIU2IP 





"11 


"12 


2 

+ 


"11 


"12 




"21 


"22 




"31 


"32 



"11 "12 


2 

+ 


"21 "22 


2 

+ 


"21 


"22 


2 
+ 


"31 


"32 


"41 "42 




"31 "32 




"41 


"42 




"41 


"42 



where || || is used also for denoting the usual Euclidian norm of a vector. 
In a similar way we conclude that the norm of the 3-vector produced by 

Ul, U2, U3 



Ui A U2 A U3 



"11 "12 "13 
"21 "22 "23 
"31 "32 "33 
"11 "12 "13 
"31 "32 "33 ys + 
"41 "42 "43 



"11 "12 "13 
yi + U21 "22 "23 
"41 "42 "43 
"21 "22 "23 
"31 "32 "33 Yi 
"41 "42 "43 



y2+ 



IS 



|ui A U2 A U3I 



"11 "12 "13 
"21 "22 "23 
"31 "32 "33 



"11 "12 "13 
"21 "22 "23 
"41 "42 "43 



||Ui||" Ul • U2 Ul 
U2 • Ul IIU2IP U2 



U3 



Ul U3 • U2 I 

"11 "12 "13 
"31 "32 "33 
"41 "42 "43 



U3 

U3 

l|2 



"21 "22 "23 
"31 "32 "33 
"41 "42 "43 



while the norm of the 4-vector produced by Ui , U2 , U3 , U4 

Ul A U2 A U3 A U4 = |U| xi 

is given by 



Ul A U2 A U3 A U4 1 



|ui| 



Ul • U2 Ul • U3 Ul • U4 

|2 



U2 ■ Ul IIU2II' U2 • U3 U2 • U4 
2 



U3 ■ Ul U3 • U2 IIU3I 
U4 ■ Ul U4 • U2 U4 • U3 



U3 • U4 

I|U4||2 



lUI 
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